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1. INTRODUCTION 


This report covers technical progress during the second year of the contract 
entitled "Determination of Coronal Magnetic Fields from Vector Magnetograms," 
NASW-4728, between NASA and Science Applications International 
Corporation, and covers the period January 1, 1993 to December 31, 1993. Under 
this contract SAIC has conducted research into the determination of coronal 
magnetic fields from vector magnetograms, including the development and 
application of algorithms to determine force-free coronal fields above selected 
observations of active regions. The contract began on June 30, 1992 and has a 
completion date of December 31, 1994. This contract is a continuation of work 
started in a previous contract, NASW-4571, which covered the period November 
15, 1990 to December 14, 1991. 

During this second year we have concentrated on studying additional 
active regions and in using the estimated coronal magnetic fields to compare to 
coronal features inferred from observations. The following is a summary of our 
research during the past year: — 

• We have redone the force-free field (FFF) fit to active region 5747 on 
October 20, 1989 using the (new) non-periodic version of the code; 

• We have documented the coronal field structure of AR5747 and 
compared it to the morphology of footpoint emission in a flare, 
showing that the "high-pressure" Ha footpoints are connected by 
coronal field lines (Mikic & McClymont 1994); 

• We have shown that the variation of magnetic field strength along 
current-carrying field lines is significantly different from the 
variation in a potential field, and that the resulting near-constant 
area of elementary flux tubes is consistent with observations of 
coronal loop in soft X-rays (McClymont & Mikic 1994); 

• We applied the evolutionary technique to an exact analytic force-free 
field of Low and Lou (1990), showing that the scheme converges 
when data for an exact field is supplied (Mikic & Barnes 1994); 

• We generated realistic models of force-free fields which model active 
regions. These models will be used to study the performance of the 
evolutionary technique to measurement errors and inconsistencies 
in the data, as well as to study the effect of photospheric shear on 
active regions, and their possible role as flare triggers (Mikic 1993). 

• We deduced coronal fields for two additional vector magnetograms, 
of AR6919 on 15 November, 1991, and of AR7260 on 18 August, 1992. 
These regions had coincident soft X-ray Yohkoh observations. 

• In conjunction with an NSF-funded grant, we studied the theoretical 
properties of an arcade field in spherical geometry. We found that 
photospheric shear can lead to magnetic nonequilibrium, causing an 
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arcade to disrupt. This effect may explain how coronal mass ejections 
are initiated (Mikic & Linker 1994). 

During this year we have presented our research results at the following 
scientific meetings, workshops, and seminars: Seminar at the University of 
Hawaii, April 1994; Invited talk at the Physics Computing '93 Conference, 
Albuquerque, June, 1993; Talk and Poster at the AAS/Solar Physics Division, Palo 
Alto, July, 1993; Invited talk at the Gordon Conference, Plymouth, New 
Hampshire, August 1993; Talk at the Sacramento Peak Workshop on Active 
Region Evolution, Sunspot, August /September, 1993; Invited talk at the 
American Geophysical Union Conference, San Francisco, December 1993. 

All the vector magnetograms used so far have been obtained at Mees Solar 
Observatory (MSO) of the University of Hawaii using the Haleakala Stokes 
Polarimeter. This project has benefited tremendously from cooperation with the 
staff of the Institute for Astronomy, University of Hawaii. Dr. Alexander 
McClymont has been our point-of-contact and principal collaborator during this 
project. 

Section 2 of this report contains an brief account of progress during the 
research performed under this contract. Section 3 contains the proposed 
statement of work for the third year of the project. 
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2. TECHNICAL PROGRESS 


2.1 Active region 5747 on October 20, 1989 

As described in last year's progress report, we have developed a non- 
periodic version of the computer code which implements the evolutionary 
technique. This represents a major advance in the realism with which we can 
model the coronal magnetic field. Previously, for numerical convenience, we 
had assumed that the transverse dimensions (the x-y plane, which is parallel to 
the photosphere) were periodic. We have been able to remove this constraint by 
using a more sophisticated numerical algorithm. We have redone our estimate 
of the coronal field for AR5747 using this new code. The analysis described below 
is based on this improved estimate. 

2.2 The Coronal Field above AR5747: Comparison with Observations 

We have made use of the coronal field computed from the magnetogram 
of AR5747 in two studies of relevance to flare and active-region physics. First, we 
compare the connectivity implied by the coronal field with Ha data obtained at 
Mees Solar Observatory (University of Hawaii) during an M2 flare (Leka et al. 
1993) which occurred 3 hours after the magnetogram scan was completed. 
Secondly, we studied the area variation along current-carrying elementary flux 
tubes, and showed that the presence of current in the corona can account for the 
observed constant thickness of soft X-ray loops (McClymont & Mikic 1994). 

Figure 1 shows observed Ha flare features superimposed on the 
magnetogram of AR5747. The three small pseudo-circular features in Fig. 1(a) 
enclose areas in which specific Ha signatures were detected. The small feature on 
the neutral line at position (3.0,5.0) showed the signature of particle precipitation 
and coincided temporally with hard X-ray emission. Figs. 1(b) and ( c ) (refer to the 
caption) show that the footpoints of a coronal loop in the estimated force-free 
coronal field are in much closer agreement with the observed high-pressure 
regions than footpoints of a loop in the potential magnetic field (Miki6 & 
McClymont 1994). 

It is unfortunate that there are no soft X-ray observations of AR5747 to 
confirm the existence of a high-pressure coronal loop joining the Ha footpoints. 
However, we have made use of this dataset in a generic way, by examining the 
variation of cross-sectional area along elementary flux tubes which carry 
substantial currents through the corona [principally from the positive spot at (3,4) 
to the negative spot at (6,8) in Figure 1(a)], and comparing the results with recent 
soft X-ray observations. Since the Skylab era it has been noted that images of 
coronal loops seen in soft X-rays or extreme ultraviolet emission appear to be of 
remarkably uniform thickness, a result recently confirmed by high-resolution 
observations (Golub 1991; Klimchuk et al. 1992). If the emitting volume outlines 
a bundle of field lines, or flux tube, the implication is that the magnetic field does 
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Figure 1. — ( a ) Ha flare features superimposed on the magnetogram of AR5747 
(from Leka et al. 1993). The neutral line is marked by a broken line; regions of 
strong vertical current are shown as shaded contours, and the three small 
pseudo-circular features enclose areas in which specific Ha signatures were 
detected. The features marked "A" at position (2.9,5.6) and "B" at (2.2,3.4) 
mark sites of high pressure "footpoints." ( b ) Field line traces in the estimated 
force-free coronal field, showing that the high-pressure sites A and B are close 
to the footpoints of a coronal loop. ( c ) Field line traces in the potential coronal 
field, with the same initial footpoint positions (at A) as in ( b ), showing that the 
agreement between the high-pressure sites and coronal loop footpoints is 
much better for the estimated force-free field than for a potential field. The 
contours in (b) and (c) show the vertical current density in the photosphere. 
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not expand in height as it would if the field were potential. Recently Klimchuk et 
al. (1992) have quantified this characteristic, finding for several loops observed by 
the Soft X-ray Telescope on the Yohkoh spacecraft a thickness variation along 
their lengths of only 10-20%. We have demonstrated that this observation is 
consistent with the characteristics of current-carrying field lines in a highly 
sheared active region. By tracing field lines in the computed coronal field for 
AR5747 we have shown that magnetic loops which are highly sheared do not 
expand rapidly in height, as they would in a potential field (McClymont & Mikic 
1994; this preprint is attached in Appendix A). 

2.3 Low and Lou's Exact Force-Free Field 

We have verified that the evolutionary technique can correctly estimate 
the coronal field when the boundary data corresponding to an analytic force-free 
field is used. The boundary data from the analytic force-free field of Low and Lou 
(1990) was used to specify B z and J z at z = 0 in the evolutionary technique, which 
was then used to estimate the coronal field. Figure 2 shows a comparison 
between the field lines in the analytic solution and those in the coronal field 
determined from the boundary data. This simulates the application of the 
evolutionary technique to a vector magnetogram, albeit to one with perfectly 
consistent data. The good agreement between the estimated coronal magnetic 
field and the analytic solution verifies that the evolutionary technique is based on 
a well-posed formulation (Mikic & Barnes 1994). 

2.4 Realistic Force-Free Field Models of Active Regions 

In order to study the effect of inconsistencies in the boundary data on the 
coronal field estimates obtained with the evolutionary algorithm, it is first 
necessary to create exact force-free coronal fields against which the estimates can 
be compared. We require models which have the properties of observed active 
region fields (complexity, no symmetries, compact flux and current distributions, 
three-dimensional variation). We have created an idealized model of a force-free 
field whose properties are similar to those of active region 5747 of 20 October 1989 
(Mikic 1993). This solution will be used for the additional purpose of constructing 
a realistic model of an active region to study flare triggers. 

We first construct a flux distribution B z (x,y ) in the photosphere which 
captures the essential features seen in the vector magnetogram of AR5747. We 
construct the vertical flux B z by superimposing four Gaussian distributions to 
match the two sunspots with negative flux and the one sunspot with positive 
flux; the net flux is zero. This flux distribution is used to compute a potential 
field in the corona. The flux distribution and the corresponding potential field 
lines are shown in Figure 3(a). 

A force-free field was generated by subjecting this potential field to specified 
footpoint displacements. The dynamical MHD equations were solved using our 
3D resistive MHD code. An equilibrium was found by introducing finite 
viscosity, a small amount of resistivity, and relaxing the field after the footpoints 
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Low and Lou Analytic 
Force-Free Field 

(a) Field Lines in Exact Field 



(b) Field Lines in Estimated Field 



Figure 2. — Comparison between field lines in (a) the analytic force free field of 
Low & Lou (1990), and ( b ) the coronal field estimated using the evolutionary 
technique from the boundary data at z = 0. The contours show the vertical 
magnetic field. The good agreement shows that the evolutionary technique is 
based on a well-posed mathematical formulation. 
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of the field lines were displaced. The displacement was introduced by applying 
flow on the photospheric boundary (which was ramped on and off gradually). 
The choice of footpoint displacements (i.e., the applied photospheric flow profile) 
was guided by the observed topology of the estimated coronal force-free field. We 
chose a flow profile to twist the fields in a manner which reproduced the 
observed global twist in the coronal field. Figure 3(b) shows the applied flow 
profile in the photosphere. 

The resulting field relaxes to a force-free field with an energy that is 40% 
above that of the potential field, consistent with that estimated for AR5747. 
Figure 4 shows a comparison between the observed vector magnetogram data and 
the photospheric profiles from this model. The vertical current density profile J z 
is remarkably similar to that deduced from the observations. In addition, the 
photospheric values of the transverse magnetic fields B x and By are similar to 
those from the observations. Figure 5 shows the field lines for the model field 
and those for the estimated coronal field for AR5747. It is apparent that the global 
features seen in the model are similar to those observed in the data. This is 
therefore quite a realistic model of the field above an active region. 

2.5 Active Region 6919 on 15 November, 1991 

We have determined the coronal field above NOAA active region 6919 on 
15 November, 1991, a region of intense interest with an extensive database of 
coincident Yohkoh and MSO observations. Figure 6 shows a projection of traces 
of the estimated coronal field above this active region, compared to the potential 
field. Note the complicated nature of this active region, which consists of several 
sunspots. The topology of the computed coronal field is consistent with the 
morphology deduced from soft and hard X-ray Yohkoh images during an X class 
flare at ~ 22:40 UT (McClymont & Mikic 1993) 

2.6 Active Region 7260 on 18 August, 1992 

We have determined the coronal field above NOAA active region 7260 on 
18 August, 1992. This region was observed with Yohkoh as it made its transit 
across the solar disk during the period August 9-22. It was associated with several 
flares, and had evidence of emerging flux. This active region was the subject of 
intense coordinated study during a recent workshop in Hawaii (November 29- 
December 3, 1993). Our calculations of the coronal field were presented at this 
workshop. 

Figure 7 shows a projection of traces of the estimated coronal field above 
this active region, compared to the potential field. Note that the field lines are 
significantly nonpotential. 

2.7 Disruption of Coronal Magnetic Arcades 

The ideal and resistive properties of isolated large-scale coronal magnetic 
arcades were studied using axisymmetric solutions of the time-dependent 
magnetohydrodynamic (MHD) equations in spherical geometry. We examined 
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(a) Potential Magnetic Field Lines 



(b) Applied Photospheric Flow 
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i- 


Figure 3. — (a) Contours of the vertical magnetic field B z {x,y) in the 
photosphere 2 = 0, with traces of the potential magnetic field lines for an 
idealized model of an active region. ( b ) Projection of the transverse 

photospheric shear velocity that is applied to the magnetic field footpoints at 
z = 0 in (a) to create a force-free field. 
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A Model Force-Free Field 
Matching Magnetic Fields Measured 
with a Vector Magnetograph 

(a) Bz, AR5747 (20 Oct.1989) Jz, AR5747 (20 Oct.1989) 

Vector Magnetograin Vector Magnetograin 



(b) Bz, Model Force-Free Field Jz, Model Force-Free Field 








min 0 max 


Figure 4. — Comparison of the normal magnetic field and current density at 
z = 0 between (a) the vector magnetogram of AR5747 on 20 October 1989 and 
(b) an idealized force-free field model. The similarity between the main 
features in the magnetogram and the model indicate that this is a realistic 
model of an active region. 
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(a) Estimated Coronal Field, AR5747 



(b) Model Force-Free Field 



Figure 5. — Traces of the field lines in (a) the estimated coronal field of active 
region 5747, and (b) the idealized force-free model of this active region. Note 
that the model field is qualitatively similar to the field of AR5747. 
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AR6919, 15 November, 1991 

(a) Potential Field 
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Figure 6. — Traces of the field lines in (a) the potential coronal field of active 
region 6919 on 15 November, 1991, and ( b ) the estimated force-free field as 
determined from a vector magnetogram of this region. 
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AR7260, 18 August, 1992 


(a) Potential Field 



(b) Force-Free Field 



Figure 7. — Traces of the field lines in ( a ) the potential coronal field of active 
region 7260 on 18 August, 1992, and ( b ) the estimated force-free field as 
determined from a vector magnetogram of this region. 
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how flares and coronal mass ejections may be initiated by sudden disruptions of 
the magnetic field. The evolution of coronal arcades in response to applied 
shearing photospheric flows indicates that disruptive behavior can occur beyond 
a critical shear. The disruption can be traced to ideal MHD magnetic 
nonequilibrium. The magnetic field expands outward in a process that opens the 
field lines and produces a tangential discontinuity in the magnetic field. In the 
presence of plasma resistivity, the resulting current sheet is the site of rapid 
reconnection, leading to an impulsive release of magnetic energy, fast flows, and 
the ejection of a plasmoid. These results are related to previous studies of force- 
free fields and to the properties of the "open-field" configuration. The field lines 
in an arcade are forced open when the magnetic energy approaches (but is still 
below) the open-field energy, creating a partially open field in which most of the 
field lines extend away from the solar surface. Preliminary application of this 
model to helmet streamers indicates that it is relevant to the initiation of coronal 
mass ejections. Appendix B contains a preprint of this article (Mikid & Linker 
1994 ). 
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3. STATEMENT OF WORK (3RD YEAR) 


During the third year of this project we intend to continue development 
and application of the "evolutionary technique" for the determination of coronal 
magnetic fields. Next year's plan is centered around the following principal goals: 

3.1 Comparison of Estimated Fields with Coronal Observations 

We will continue to compare the properties of the estimated coronal 
magnetic fields with observations. Initially we will compare the magnetic fields 
already computed for AR6919 on 15 November, 1991, and AR7260 on 18 August, 
1992 (as described in Sections 2.5 and 2.6) to Yohkoh soft X-ray images of these 
regions, and also to Yohkok hard X-ray emissions and Ha brightening during 
flares. 

3.2 Application to Additional Vector Magnetograms 

We are in the process of using the excellent data from the Advanced Stokes 
Polarimeter (ASP), which was built at HAO, and is presently used for solar 
observations at the Sacramento Peak Observatory in Sunspot, New Mexico. This 
instrument has very high spatial resolution, and is a "next-generation" 
polarimeter with the capability of measuring photospheric and chromospheric 
magnetic fields with high accuracy. We will be analyzing magnetograms taken in 
June 1992, and also those taken recently in November 1993. 

3.3 Effect of Emerging Flux 

We will continue a theoretical effort, which was begun this Fall, to study 
the effect of emerging flux. Significant emergence of magnetic flux was observed 
during the evolution of AR7260 in August 1992. (Recall that we have estimated 
the coronal field above this region on August 18, as described in Section 2.6.) The 
theoretical role of emerging flux is not well know. We are developing our code 
to include the effect of emerging flux, and we will study its effect during the next 
year. 

3.4 Accuracy of Computed Coronal Fields 

We will continue to study the accuracy and limitations of the evolutionary 
technique. The effect of errors will be studied by applying the technique to 
"synthesized magnetograms," which will be constructed by adding random errors 
to fields with known solutions. The solutions we will use are those described in 
Sections 2.3 and 2.4. 
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ABSTRACT 


It has been noted for many years that images of active region coronal loops seen in soft 
X-rays or extreme ultraviolet emission suggest a pipe-like appearance. Recently Klimchuk et 
al. (1992) have quantified this characteristic, finding for several loops observed by the Soft X- 
ray Telescope on the Yohkoh spacecraft a thickness variation along their lengths of only 10— 
20%. We demonstrate here that this observation is consistent with the characteristics of 
current-carrying field lines in a highly sheared active region. Vector magnetogram data on 
NOAA active region 5747, taken with the Stokes Polarimeter at Mees Solar Observatory on 
1989 October 20, provided photospheric boundary conditions from which a force-free coronal 
magnetic field was computed. By tracing field lines, we show that magnetic loops which are 
highly sheared do not expand rapidly in height, as they would in a potential field. In addition, 
the expanding sections close to the footpoints of current-carrying twisted loops tend to be more 
vertical than in a potential field, so that when seen projected against the solar disk, the loops 
appear to terminate more abruptly. Consequently current-carrying loops exhibit a near-uniform 
cross-section with thickness variations of order 30% along their lengths. 

Subject headings: Sun: Corona — Sun: Magnetic Fields — Sun: X-Rays, Gamma Rays 

1. INTRODUCTION 

Coronal soft X-ray emitting loop structures are generally believed to identify the paths of % 
magnetic field lines. Since the Skylab era it has been noted that some coronal loops appear to if 
be of remarkably uniform thickness, a result recently confirmed by high-resolution observations * 
(Golub 1991, Klimchuk et al. 1992.) If the emitting volume outlines a bundle of field lines, or ' 
flux tube, the implication is that the magnetic field does not expand in height as it would if the 
field were potential. An alternative interpretation is that the emitting volume does not outline a 
flux tube, but that emission near the footpoints of the loop encompasses a thicker bundle of 
field lines than emission from the upper part of the loop. The temperature/density structure of 
a loop in energetic equilibrium (Craig, McClymont & Underwood 1978, Rosner, Tucker & 
Vaiana 1978) could reduce the diameter of the emitting region which is above the threshold for 
detection as height increases. However, it seems very unlikely that the distribution of emis- 
sivity along the loop would compensate so precisely for the change in magnetic cross-section, 
or that the apparent constant thickness of the loop would be preserved in observations sensitive 
to different temperatures. Here we assume that the emission outlines a flux tube, and show 
that the coronal field of a highly sheared active region, reconstructed from actual vector mag- 
netograph data, is consistent with this interpretation. 

The variation of magnetic field strength along a flux tube is of interest from several other 
points of view. One is the effect of area variation on the thermal conductive energy balance of 
the plasma in a quasi-static hot loop (e.g.. Dowdy, Emslie & Moore 1987). Area variation 
may also have important effects on mass flows, particularly in the context of blue-shifted soft 
X-ray emission in the impulsive phase of flares, which can be interpreted as the signature of 
chromospheric evaporation (see Doschek 1990, Antonucci, Dodero & Martin 1990). Another 
important consequence of area variation is the magnetic mirroring of high energy particles 
accelerated in the impulsive phase of a flare (e.g. Petrosian 1985, LaRosa & Emslie 1988, Lu 
& Petrosian 1990). Theoretical investigations of these effects to date have adopted simple ad 
hoc models of the area variation. The results of these studies also seem to be contradictory to 
some extent. For instance, LaRosa & Emslie (1988) conclude that to explain the low bright- 
ness of the chromosphere relative to coronal emission in the impulsive phase of a flare, a loop 
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into which non thermal electrons are injected must narrow considerably towards its footpoints. 
On the other hand, Petrosian (1985) suggests that the field must be nearly uniform to explain 
impulsive y-ray directivity. 

We present here the first theoretical study of area variation along magnetic flux tubes 
based on observational data. The reduction of the vector magnetograph data is described in § 
2, and the interpretation of these data in terms of photospheric boundary conditions on the 
coronal magnetic field is discussed in § 3. The method for constructing a force-free represen- 
tation of the coronal field is described in § 4, and the analysis of loop thickness variations is 
presented in § 5. Our results are described in § 6, and our conclusions summarized in § 7. 

2. OBSERVATIONS 

The data used in this analysis were obtained with the Stokes Polarimeter at Mees Solar 
Observatory on Haleakala (Mickey 1985). NOAA active region 5747, which produced a series 
of major flares, was observed 1989 October 18 — 22. Here we use the data obtained on 
October 20 (scanned from 17:41 to 18:33 UT), which are of superior quality compared to the 
other days. 

A full account of the observations and data reduction techniques has been given by 
Canfield et al. (1993) and Leka et al. (1993). Briefly, the Stokes Polarimeter obtained full 
Stokes profiles in the magnetically sensitive Fe I lines at 6301.5 and 6302.5 A, building up a 
raster by scanning poim-by-point over the active region, in steps of 5.6", with a field of view 
at each pixel of 6". The data were reduced using the least-squares fitting code of Skumanich 
& Lites (1987). Raster points at which the polarization was too low for a least-squares fit 
were reduced using the integral method (Ronan, Mickey & Orrall 1987). The ambiguity in 
direction of the transverse magnetic field in the resulting vector magnetogram was resolved as 
described by Canfield et al. (1993), and the magneto gram transformed to disk center. The 
vertical magnetic field and current density, computed from the curl of the horizontal field, 
define the photospheric boundary conditions needed to reconstruct the coronal force-free field. 

We have also used observations of a near-potential active region for comparison with the 
highly nonpotential AR 5747. These observations, of AR 7490, magnetic classification a, 
were scanned with the Haleakala Stokes Polarimeter on 1993 April 28 from 17:31 to 19:03 
UT, and reduced using the same procedure outlined above. 

3. INTERPRETATION OF OBSERVATIONS 

Since the photosphere is probably not force-free at the observed level (i.e. the depth of 
formation of the spectrum lines used in the analysis), is it then possible to use this data to 
compute a force-free field? The effect of mechanical forces can be considered in the limit of 
two spatial scales: (1) individual flux tubes which are much smaller than both the spatial reso- 
lution of the instrument and the global size scales of the active region, and (2) the large scale 
field of sunspots. 

Outside sunspots, the magnetic field at the photosphere does not vary smoothly, but 
occurs in flux tubes of high field strength (=1500 G) but very small size scale (=150 km), 
separated by field-free photosphere (Stenflo 1989). Therefore the field lines fan out rapidly as 
they leave the pressure-confining photosphere. Since the field spreads out to become force-free 
only a few hundred km above the photosphere, the discrete nature of the field at the lower 
level should not affect the global structure of the corona. The “spread-out” field is the 
appropriate one to use as a boundary condition in computing the coronal field. The question is 
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whether we can approximate this field using the data available. Since the spatial resolution of 
the polarimeter, 4-5 Mm, is much coarser than the fine scale photospheric structure and is of 
order the smaller size scales important to the global structure of the active region, the polarime- 
ter does indeed yield a suitably smoothed field, provided it correctly measures the average field 
strength. The circular polarization signal is proportions' ^ : longitudinal field strength, and 

so measures the average field directly. The linear polarization signal, however, is proportional 
to the square of the transverse field, and thus overestimates the average transverse field when 
the filling factor is low (Stenflo 1985). In strong field regions (sunspot umbrae), this is not a 
problem, since the filling factor is very close to unity. Over most of the area of significant flux 
in the active region (penumbrae, plage), we can use the information in the line profiles to 
determine the filling factor (Skumanich & Lites 1987), so that the average transverse field over 
a pixel can be derived, at least to first order. In weak field regions (where there is little visual 
indication of the presence of magnetic field), the polarization signal is too small for a least- 
squares fit, and we use the integral method (Ronan, Mickey & O trail 1987), which does not 
yield a correction for the filling factor. Since the filling factor is generally small in these areas, 
the transverse field may be overestimated. However, the contribution to the active region flux 
from these weak field regions is of order 1% in the present dataset, so these errors will not 
have an appreciable effect on the accuracy of tire force-free solution. 

Turning to the large scale mechanical forces present in the photosphere, i.e. the gas pres- 
sure confining sunspots, the distribution of the photospheric magnetic field must satisfy a 
number of constraints to be compatible with an overlying force-free field (e.g. Low 1985, Aly ^ 
1989). First, in order for the region to be considered isolated, the magnetic flux and vertical^ 
current must balance within the field of view. In the present case, both balance to within al 
few percent Second, the total force and torque exerted on the photosphere by the overlying * 
magnetic field must vanish. For this dataset, the horizontal forces and the vertical torque are 
of order a few percent (normalized by the mean magnetic pressure integrated over the photo- 
sphere, and in the case of torque multiplied by the size of the active region), but the vertical 
force and horizontal torque are of order 10%. These particular discrepancies are not unex- 
pected however, since the photosphere is manifestly not force-free, in the sense that the mag- 
netic field in sunspots would be expelled from the photosphere by the buoyancy force if it was 
not attached to a subphotospheric flux tube. The total vertical force on the photosphere is 
given by F n = j(B?-B, 2 ) dA /8x, where B n and B, are the magnetic field components normal to 
and parallel to the photosphere. In a sunspot, the field lines are pushed together into a more 
vertical orientation than if they were force-free, yielding a positive F „ . However, the currents 
associated with this force circulate horizontally around the flux tubes, extending at most a few 
pressure scale heights (a few hundred km) above the observed level. They do not enter die 
corona, and therefore do not influence in a direct way the global structure of the coronal mag- 
netic field. These currents cannot be detected in a vector magnetogram made at one level in 
the atmosphere^ which can measure only the vertical component (/„) of the current The non- 
force-free nature of sunspots is not expected to change the structure of the coronal field in any 
substantial way (Aly 1989). 

Another constraint which the photospheric field must satisfy if the coronal field is to be 
force-free concerns the distribution of a=/„/B„. Since a is conserved along each field line, 
there must be the same number of footpoints of field lines with a given value of a in each 
magnetic polarity. A study of this same magnetogram (Canfield n al. 199!) revealed 
discrepancies in the distributions of flux over a which seemed to be significantly above the 
noise. One of our long-term aims is to investigate the response of the numerical model to such 
discrepancies in the data. 
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4. COMPUTATION OF THE CORONAL MAGNETIC FIELD 

To compute the coronal magnetic field, a 3D resistive MHD code (Mikic, Baines & 
Schnack 1988, Mikic & Barnes 1993) is used to evolve the field until it satisfies the desired 
boundary conditions at the photosphere. The magnetogram, which does not map onto a rec- 
tangular region of the solar surface, is first interpolated onto a Cartesian grid. This is embed- 
ded in a larger field-free region to isolate the active region from the effect of the computational 
boundaries. The grid is nonuniform in the horizontal dimensions, to remove the boundaries as 
far as possible while maintaining resolution in the region of interest The grid is also nonuni- 
form in the vertical direction, to resolve rapid variations in the field close to the photosphere. 
The 3D grid (60x60x40 mesh points) forms a cuboid 366 Mm on a side, roughly four times 
the extent of the strong-field region of the spot group, and 244 Mm high. The potential field 
defined by the normal component of the photospheric field, B n , is computed first, then voltages 
are applied at the photosphere to drive currents through the corona. The voltages are ramped 
up over time and adjusted adaptively to drive the currents towards the desired normal current 
density, J n . Resistivity allows the evolving field lines to change their topology as required to 
reach a force-free state. Once the field is close to its final configuration, the resistivity can be 
reduced, and if the boundary data are perfectly consistent, a force-free field will ensue. 

In practice, data imperfections, the effect of discretization, and possibly mechanical forces 
in the photosphere, prevent a completely force-free state from being attained. We have found 
that a small remnant JxB force is present near the photosphere, although it decreases rapidly 
with height. The corresponding steady state solution has finite flows, with steady electric fields 
at the boundary to balance the resistive dissipation; the JxB force is balanced by (artificially 
large) viscous forces. In order to obtain a completely force-free field we have also tried a 
different prescription. Starting from the abovementioned steady state solution, we reduce the 
resistivity to zero while fixing the positions of the magnetic field footpoints, without applying 
boundary potentials, and allow the field to relax. We find that the field then relaxes to a 
force-free state in which the normal electric current density at the photosphere does not match 
exactly that deduced from the vector magnetogram. For the data reported here, the normal 
electric current density is within 20-30% of the magnetogram current density. The two mag- 
netic fields produced by this technique can be considered as two (equally good) approximations 
to the true field. At heights greater than =5 Mm above the photosphere, the two coronal fields 
do not differ significantly. The latter field is used in the analysis described in this paper. 

5. ANALYSIS OF LOOP THICKNESS VARIATIONS 

To extract the information used in this study from the computed coronal field, we trace 
field lines by integrating the equation drlds =B/lfl I, where r is the position vector and s is 
arclength along the field line, from selected photospheric footpoints. An adaptive-step second 
order Runge-Kutta algorithm, with trilinear interpolation to obtain the magnetic field between 
mesh points, proves to be of adequate accuracy (r.m.s. error in footpoint position of approxi- 
mately 0.03 mesh points when tested on a dipole field). 

For the purposes of this study, we adopt the simple “microscopic” definition of flux tube 
cross-sectional area, A where <b is the (infinitesimal) magnetic flux in an elementary 

flux tube. The variation in diameter (width) along an elementary flux tube, assuming, to a first 
approximation, a circular cross-section, is then given by w «* The variation of w along a 

field line (of interest for particle mirroring studies) can be studied by plotting w against s. 
The apparent thickness of soft X-ray emitting loops, seen in projection against the solar disk, 
can be studied by plotting w as a function of distance along the loop projected onto the 

•* •- '■ 


- 5 - 



photosphere. We have not evaluated quantitatively the actual width of a “macroscopic” flux 
tube. A more detailed investigation would have to act unt for the changing shape of the bun- 
dle of field lines along the length of the loop. We do. .xjwever, compare qualitatively macros- 
copic flux tubes in the highly nonpotential AR 5747 with similar macroscopic flux tubes in the 
near-potential region AR 7490. 

In the complicated three-dimensional structure of the active region magnetic field, many 
of the variables are correlated (e.g. the magnetic field falls with increasing height, strong 
currents only occur in regions of strong field, etc.), and it is difficult to demonstrate directly 
that field-aligned currents are responsible for changes in the shape of flux tubes. Here we 
adopt three approaches: 

First, we compare the thickness variation along field lines traced in the computed force- 
free field with field lines traced in the potential field computed from the same boundary data 
but ignoring currents entering the corona. We choose footpoints in the areas of strong photos- 
pheric current density. Here we deal with the morphology of only one active region, but a 
drawback of this approach is that the potential field can be criticized as unphysical: for 
instance, force-free field lines starting from footpoints in a sunspot may end in a sunspot of 
opposite polarity, while potential field lines starting from the same footpoints may end in a 
weak field region. 

Second, we have studied the correlation between current density at the footpoints of the 
field lines and the “expansion factor” (following the definition of Klimchuk et al. 1992, but 
applied to “microscopic” flux tubes). Klimchuk et al. "straighten” observed loops, seen pro- 4 
jected against the photosphere, and define the expansion factor as the ratio of apparent loop T 
thickness at the midpoint between the footpoints to the apparent thickness at the narrower foot- 
point. They define the “footpoints” to be the locations at which the soft X-ray intensity drops 
to half of the maximum intensity along the loop. Since we do not have a measure of the X-ray 
brightness of our model loops, we take the “footpoints” to be 15% of the loop length (pro- 
jected) in from the true photospheric footpoints. This corresponds to the average distan ce 
between the half-intensity points and the apparent ends of the observed loops (Klimchuk, 
private communication). Because of the correlation of strong current with strong field, and 
stronger fields with lower heights, careful interpretation of these scatter {riots is necessary. 

Third, to achieve a comparison with physically realistic potential field lines, and also to 
examine the characteristics of “macroscopic” flux tubes, we compare the shapes of flux tubes 
in the force-free field of AR 5747 (again rooted in strong current regions) with those in a 
different active region which is close to potential (AR 7490, observed with the Haleakala 
Stokes Polarimeter 1993 April 28). Although the potential field lines are physically realistic in 
this case, the active region has a very different character from AR 5747. 

6. RESULTS 

Figure 1 shows the magnetogram of AR 5747 with field lines projected against the photo- 
sphere. Contoura show the vertical component of the photospheric magnetic field. In the left 
panel, the field lines were computed using the reconstructed force-free field, while the field 
lines in the right panel were computed from the potential field obtained by ignoring the 
currents entering the corona. The same footpoint positions, chosen to lie in the areas of strong 
vertical current (in both magnetic polarities), were used in -adng boih iats of field lines. It 
will be noted that in the potential field calculation, field lines originating in the left part of the 
positive spot cross the neutral line in a perpendicular direction, while in the force-free calcula- 
tion, the same field lines are strongly sheared, running a long distance parallel to the neutral 
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line, to end in the upper of the two negative spots. 

Figure 2 shows the width (diameter) of elementary flux tubes (assuming that the cross- 
section remains circular) as a function of distance along the loop, projected onto the photo- 
sphere. The widths are normalized to unity at the footpoint with the stronger field. This simu- 
lates the measured width of soft X-ray loops seen in projection against the solar disk. On the 
left are field lines from the computed force- free field; on the right, for comparison, are field 
lines from the potential field, computed by ignoring the currents entering the corona from the 
photosphere. The difference is striking. Most of the force-free flux tubes are of remarkably 
uniform diameter over most of their lengths, narrowing only close to their photospheric foot- 
points (cf. Golub 1991). The potential field lines, in contrast, show no such attribute. The 
variation in thickness of the force-free flux tubes may be estimated by eye as being of order 
10 — 20% over most of their length, as observed by Klimchuk et al. (1992). A comparison of 
the computed shapes with the observed shape shown in Figure 2 of Klimchuk et al. (1992) 
seems compelling. 

To gain insight into the difference between the force-free and potential fields, it is useful 
to consider the variation in area along each field line as a function of arclength from each foot- 
point. Figure 3 shows such a plot, in which each curve represents the portion of a field line 
from its footpoint to its widest point. In both the force-free and potential fields, the initial rate 
of increase in diameter with distance from the footpoint is about the same up to a distance of 
=10 Mm, which corresponds roughly to the radius of the large spots. Therefore, close to the 
footpoints the decrease in field strength is controlled largely by the fanning out of field lines as 
they leave the sunspots, independent of the character of the active region on the global scale. 
At greater heights, there is a significant difference, with the diameter of current-carrying flux 
tubes increasing much less rapidly than in the potential field case. 

We have attempted to quantify the effect of field-aligned current on the shape of flux 
tubes by using the “expansion factor” as defined by Klimchuk et al. (1992). As mentioned in 
§ 5, it is difficult to unambiguously relate the shape of the flux tube to current density, since, 
for example, the strongest currents occur only in the strongest magnetic fields; because the sun- 
spot group is relatively compact in this case, these field lines are low-lying. Therefore scatter 
plots of expansion factors of selected field lines against current density at the photosphere (Fig. 
4a), and of expansion factor against maximum height attained by the field line (Fig. 4b), both 
show a correlation. The plot of expansion factor against field strength at the photosphere (Fig. 
4c) shows less correlation. In these plots, we have averaged the current densities and field 
strengths between the two footpoints, and normalized them to the peak values of approximately 
27.2 mA m~ 2 and 2840 G present in the magnetogram. The maximum heights were normal- 
ized to 50 Mm, the distance between the centers of the two largest sunspots. The set of field 
lines used here is larger than the set shown in Figure 1, and includes field lines from regions 
of low photospheric current density. We excluded from the set a few field lines with apices 
higher than the maximum height shown in Figure 4b; some of these have much larger expan- 
sion factors than shown here, and the highest tend to be distorted by the upper computational 
boundary, 

Can we conclude that current-carrying field lines are of more constant cross-section than 
non-current-carrying field lines, or should we conclude that lower-lying field lines are of more 
constant cross-section? It is clear from Figure 4a that on average the expansion factor 
decreases with increasing current density. For field lines canying significant currents (> 1/4 of 
the maximum current density), the expansion factor is below 1.5, regardless of the height the 
field line reaches. The average expansion factor in the high-current range is approximately 1.3, 
somewhat larger than the value of 1 . 1 3±0. 10 found by Klimchuk et al., but nevertheless 
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suggesting a significant result However, some loops carrying only small currents also exhibit 
small expansion factors (lower left area of Fig. 4a.) 

In Figure 4b we see that the smallest expansion factors of all occur for very lowlying 
loops ( H < 0.1, or 5 Mm.) These are short loops (average length around 15 Mm) bridging the 
neutral line. Because of their small size compared to the global scale of the active region, the 
magnetic field strength does not vary much along these loops. Their small expansion factors 
are due to their situation in a two dimensional arcade field, while higher field lines expand in 
three dimensions. Although these field lines are strongly sheared in this case, they would have 
relatively constant cross-sections even in a potential field. If the potential field of an active 
region could be represented by a single buried horizontal dipole, at depth D , (not a good 
approximation in general) the field strength at the loop tops (for field lines lying in the vertical 
plane above the dipole) would fall off with height H, as (D+H ) -3 , implying an increase in 
flux tube diameter at the apex with height of (D +H ) m . Measured along an individual field 
line, the flux tube diameter would increase linearly with height, h , except close to the apex 
(h=H): 

w°c (D+h) {1 -V4[(D+h)/(D+H)f 3 r VA . (1) 

This yields the expansion factor 

/ = 1 + (5/3)// /D H<D 

/ =V2( l+HID ) H>D ( * 

Estimating a dipole depth of 50 Mm, the distance between the biggest spots, then leads us to * 
conclude that field lines (in the vertical plane) with apex heights below 10 Mm (0.2 normalized + 
units) will have expansion factors of less than 1.33 solely because their lengths are short com- - 
pared to the scale size of the global field. This is consistent with the behavior seen in Figure 
4b. For higher-lying field lines (those with 0.2 <H <0.7, corresponding to maximum heights 
of 10 Mm to 35 Mm) we attribute small expansion factors to the presence of field-aligned 
currents. For H >0.7 (35 Mm) the lower limit on the expansion factor seems to drift upwards, 
becoming greater than 1.5 for field lines reaching heights of 50 Mm. The largest expansion 
factors shown in Figure 4b are consistent with the result of equation (2) for H ~ D . 

To clarify the relationship implied by these correlations, we have constructed a scatter 
plot of the expansion factors, /, against the linear combination of the variables, 
u ~a J+b B+c H , where a, b, and c are constants to be adjusted to minimize the scatter of 
points in a least squares sense. The function fitted to the data was chosen to be a curve 
defined parametrically by cubic polynomials u(t) and / (r). The best correlation between u 
and /, shown in Figure 4d, was found for the combination of variables, u = 

J -0.18// -0.09B. Figure 4d also shows the fitted curve. This confirms the results sug- 
gested in Figure 4a-c: the magnetic field strength at the footpoint does not greatly influence the 
expansion factor; the height reached by the field line has a greater effect; but the primary deter- 
minant of expansion factor is the current density. 

So far we have considered elementary flux tubes, assumed to be of circular cross-section. 
Lastly, we compare "macroscopic” flux tubes which carry substantial currents in the highly 
nonpotential active region AR 5747, with similar flux tubes traced in a near-potential active 
region, AR 7490. The differing appearances of the “loops” shown in Figure 5 supports our 
assertion that field-aligned currents are responsible for the near-invariant cross-sections of soft 
X-ray loops. The loops shown in the nonpotential active region probably have expansion fac- 
tors of no more than 2, while the expansion factors of the loops in the potential active region 
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are cleaiiy greater than 2. 


7. CONCLUSIONS 

The intriguing observation that many soft X-ray emitting loops, whose emissions prob- 
ably outline bundles of magnetic field lines, appear to be of near-constant cross-section, can be 
explained if the loops lie on strongly sheared, current-canying field lines. Not only does the 
magnetic field strength fall less rapidly with height when current is present in the corona, but 
the twist in the field lines tends to make the force-free field lines enter the photosphere at a 
steeper angle, foreshortening the tapered section of the flux tubes seen in projection against the 
solar disk. For the dataset examined here, the variations in loop thickness are of order 30%. 
While larger than the 13% variation observed by Klimchuk et al. (1992), the variations are 
significantly less than those of a potential field. We find also that field lines in the nonpoten- 
tial active region which do not carry significant currents exhibit area variations more like the 
potential field than like the current-carrying field lines illustrated here. Thus the observation 
that bright soft X-ray loops tend to be of uniform width suggests, as mentioned by Klimchuk 
et al., that coronal heating is associated with strong field-aligned currents. Lowlying loops 
which are short compared to the global length scales of the active region form an exception to 
this statement, since they can have near-constant field strength along their lengths. A goal for 
the future will be to model the coronal magnetic field in an active region for which cotemporal 
soft X-ray observations are available. We anticipate that the techniques of modeling 3D 
force- free coronal fields employed here will answer many other questions regarding active 
region magnetic topology and energetics. 

This research was supported by NSF grant ATM-9108369 to SAIC and by NSF grant 
ATM-9106052 to tire University of Hawaii. We thank Tim Klimchuk for helpful information 
and comments, and Bob Rosner for constructive criticisms which have greatly improved the 
analysis. 
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FIGURE CAPTIONS 

Fig. 1 — Magnetogram of NO A A active region 5747 on 20 October 1989, after resolution of 
the ambiguity in the azimuth angle of the transverse field and rotation to disk center. Solar 
North is approximately 20° anticlockwise from the top. Contours (at approximately 160x2'* 
Gauss, n=0,l,2,3,4) show the vertical component of the magnetic field, with solid lines mark- 
ing positive polarity and broken lines negative polarity. In the left panel field lines with foot- 
points in areas of high photospheric current density are shown projected against the photo- 
sphere. In the right panel field lines with the same starting footpoints, traced in the potential 
field obtained by ignoring the currents entering the corona, are shown for comparison. 

Fig. 2 — Diameters of elementary flux tubes along each of the field lines shown in Figure 1, 
as a function of distance along the field lines, projected onto the photosphere. The left panel, 
showing the field lines computed from the reconstructed force-free field, illustrates the widths 
of coronal loops seen in projection against the solar disk (if the active region was at disk 
center). The right panel shows the diameters of flux tubes on the potential field lines with the 
same set of starting footpoints. While the potential field exhibits the expected variation, the 
force- free flux tubes are of remarkably constant thickness over most of their lengths (cf. Fig. 2 
of Klimchuk et al. 1992). 

Fig. 3 — Diameters of flux tubes as a function of arclength, measured along the field lines 
from each footpoint, for the force-free field (left panel) and the potential field (right panel). 
The potential field flux tubes clearly expand much more rapidly than the force-free flux tubes. 

p ig- 4 — Expansion factor (the ratio of flux tube width at the apex of the loop divided by the 
width near the narrower footpoint), for a set of field lines traced in the nonpotential force-free 
field, with a range of magnetic field strengths and current densities at their footpoints, and a 
range of apex heights. Current densities and field strengths represent averages over the two 
footpoints, and are normalized to peak values of 27.2 mA m -2 and 2840 G. The apex heights 
are normalized to 50 Mm, the distance between the main spots. The expansion factor 
decreases with increasing footpoint current density [panel (a)], and increases with increasing 
apex height [panel (b)J, but there is little correlation between expansion factor and field 
strength [panel (c)]. The combination of these parameters shown on the abscissa of panel (d) 
gives, in a simple linear model, the tightest correlation. The line represents the least-squares fit 
of a cubic parametric curve, [u(r)/(r)}. For current densities > 1/4 of the maximum foot- 
point current density in the active region, the loops expand only about 30%. While this is 
larger than the expansion factor of 13±10% found by Klimchuk et al (1992), it is considerably 
less than expected if the field were potential. Small expansion factors are also seen for lowly- 
ing field lines in panel (b), even where the current density is small [the points to the lower left 
in panel (a)]. 

F>g- 5 — “Macroscopic” magnetic loops, formed by tracing bundles of field lines in (a) the 
strongly nonpotential active region AR 5747, and (b) the near-potential active region AR 7490, 
superposed on the photospheric magnetograms showing contours of the vertical magnetic field 
strength. The nonpotential loops appear to have a more uniform thickness than the potential 
loops. 
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ABSTRACT 


The ideal and resistive properties of isolated large-scale coronal 
magnetic arcades are studied using axisymmetric solutions of the time- 
dependent magnetohydrodynamic (MHD) equations in spherical geometry. 
We examine how flares and coronal mass ejections may be initiated by 
sudden disruptions of the magnetic field. The evolution of coronal arcades in 
response to applied shearing photospheric flows indicates that disruptive 
behavior can occur beyond a critical shear. The disruption can be traced to 
ideal MHD magnetic nonequilibrium. The magnetic field expands outward 
in a process that opens the field lines and produces a tangential discontinuity 
in the magnetic field. In the presence of plasma resistivity, the resulting 
current sheet is the site of rapid reconnection, leading to an impulsive release 
of magnetic energy, fast flows, and the ejection of a plasmoid. We relate these 
results to previous studies of force-free fields and to the properties of the 
"open-field" configuration. We show that the field lines in an arcade are 
forced open when the magnetic energy approaches (but is still below) the 
open-field energy, creating a partially open field in which most of the field 
lines extend away from the solar surface. Preliminary application of this 
model to helmet streamers indicates that it is relevant to the initiation of 
coronal mass ejections. 

Subject headings: MHD — Sun: magnetic fields — Sun: corona — Sim: flares 


Suggested running head: DISRUPTION " CORONAL ARCADES 



3 


1. INTRODUCTION 

Coronal mass ejections and solar flares are spectacular manifestations 
of solar activity. It is believed that they are initiated by the sudden release of 
energy stored in the coronal magnetic field. The magnetic energy in the 
corona can be built up through shear introduced by photospheric motions. In 
principle, the free magnetic energy (that in excess of the energy in the 
potential field) is available for release. Although there is ample observational 
evidence for the existence of highly non-potential coronal magnetic fields in 
active regions ( e.g ., Gary et al. 1987; Hagyard 1988; Canfield et al. 1993; Leka et 
al. 1993), which consequently have considerable free magnetic energy, it has 
been difficult to demonstrate theoretically that this energy can be released 
impulsively. In addition, since coronal mass ejections drive solar material 
out of the corona, at least some of the magnetic field lines must also be 
opened. 

Many investigations have attempted to explain how energy can be 
released from the coronal magnetic field (Low 1981; Birn & Schindler 1981; 
Zwingmann 1987; Priest 1988; Priest & Forbes 1990; Forbes & Isenberg 1991; 
Forbes 1992; Isenberg, Forbes, & D6moulin 1993), including studies of force- 
free equilibria (Barnes & Sturrock 1972; Yang, Sturrock, & Antiochos 1986; 
Klimchuk, Sturrock, & Yang 1988; Klimchuk & Sturrock 1989; Finn & Chen 
1990; Porter, Klimchuk, & Sturrock 1992; Klimchuk & Sturrock 1992; 
Roumeliotis, Sturrock, & Antiochos 1993), the dynamical evolution of 
arcades (Mikic, Barnes, & Schnack 1988; Biskamp & Welter 1989; Forbes 1990; 
Steinolfson 1991; Wu et al. 1991; Finn, Guzdar, & Chen 1992; Inhester, Birn, & 
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Hesse 1992; De Vore & Antiochos 1993), the equilibrium properties of open 
and partially open fields (Wolf son & Low 1992; Wolfson 1993; Low & Smith 
1993), and the asymptotic properties of sheared coronal fields (Aly 1984, 1985, 
1988, 1990). However, a basic problem remains unresolved. Aly (1984, 1991) 
and Sturrock (1991) have found that the energy in an open-field 
configuration 1 is an upper limit to the magnetic energy for all force-free 
equilibria. 2 Therefore, an impulsive transition from a closed, force-free 
cc figuration to an open field would appear to be energetically impossible. 
How then can a phenomenon such as a coronal mass ejection be initiated by 
magnetic energy release? 

To answer this question we have investigated the dynamical evolution 
of an isolated arcade in axisymmetric spherical geometry using the ideal and 
resistive magnetohydrodynamic (MHD) models. Although this is an 
idealized model problem, it significantly extends previous models of coronal 
fields, and has allowed us to answer several important theoretical questions 
related to the disruption of coronal arcades. When a dipole field is subjected 
to a prescribed photospheric shear flow profile, its magnetic energy increases 
until it approaches the open-field energy, at which point the field 
configuration becomes very sensitive to additional shear, expanding 
considerably and producing a concentration of the electric current density. 
Past a critical level of shear, when the magnetic energy approaches, but is still 

1 An open field is a configuration in which one end of every field line intersects the 
photosphere, the other end extending out to infinity (e.g., Barnes Sc Sturrock 1972; Sturrock 
1991, and Section 5.2 of this paper). 

2 When the configuration has field lines which are not connected to the photosphere, the force- 
free field energy can exceed the open-field energy (e.g., Priest & Forbes 1990). 
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below, the open-field energy, the configuration experiences ideal MHD 
magnetic nonequilibrium, leading to an opening of the field and the 
formation of a tangential discontinuity in the magnetic field. The 
assumption of quasi-static behavior (in which the field is assumed to evolve 
through a series of equilibrium states) breaks down when the critical shear is 
exceeded, since the equilibration time for the expanding field becomes long, 
requiring a dynamical model to describe the system properly. The appearance 
of a tangential discontinuity implies that even a small plasma resistivity 
becomes important. The resistive MHD evolution shows that finite 
resistivity resolves the tangential discontinuity into a current sheet at which 
there is rapid reconnection of the magnetic field, leading to the release of 
magnetic energy, fast flows, and the ejection of a plasmoid, demonstrating 
that magnetic energy can be released impulsively in a coronal arcade. The 
opening of the field, and the subsequent reconnection, have been predicted 
on the basis of the asymptotic properties of sheared force-free arcades (Aly 
1985, 1990). 

There is considerable confusion in the literature in the terminology 
that has been used to describe the phenomenon of "absence of MHD 
equilibrium." The terms "loss of equilibrium" and "magnetic 
nonequilibrium" have become popular. Although they are used 
synonymously at times, we wish to distinguish between them. We refer to 
magnetic nonequilibrium as the ideal MHD process by which a plasma with 
an initially smooth magnetic field evolves into a configuration with 
inevitable tangential discontinuities in the magnetic field (Parker 1972, 1979; 
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Priest 1981; Moffatt 1985; Vainshtein & Parker 1986; Vainshtein 1987 3 ). In the 
presence of finite plasma resistivity these tangential discontinuities are 
resolved into thin layers called current sheets ( e.g ., Parker 1986). For 
simplicity we shall loosely refer to both the tangential discontinuities which 
occur in ideal MHD and the current sheets which occur in resistive MHD 
simply as current sheets. The term loss of equilibrium has been used to 
describe the disappearance of solutions to the equilibrium equations in 
response to the variation of an external parameter. [The external parameter 
could define the magnetic field or pressure profile (e.g., Low 1977; Birn & 
Schindler 1981), or the displacement of the magnetic field footpoints (e.g., 
Priest & Milne 1980).] Loss of equilibrium is characterized by the 
disappearance of a local extremum in the potential energy of the system when 
the external parameter reaches a critical value. While this "critical" behavior 
has been interpreted as an eruption of the configuration (Bim & Schindler 
1981), it does not necessarily have physical significance (Aly 1985; Klimchuk & 
Sturrock 1989; Finn & Chen 1990). 

In the definition of loss of equilibrium there is no explicit reference to 
the smoothness of the magnetic field. It is this property that distinguishes 
magnetic nonequilibrium from loss of equilibrium. Since it is necessary to 
know the potential energy surfaces of the configuration to determine whether 
loss of equilibrium has occurred, in our dynamical approach it is not possible 
to explicitly identify loss of equilibrium. On the other hand, since we can 
detect when a current sheet forms, it is possible to identify the onset of 

^ Although this paper illustrates magnetic nonequilibrium, the phenomenon has unfortunately 
been translated from the original Russian as loss of equilibrium ! 
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magnetic nonequilibrium. It is conceivable that magnetic nonequilibrium 
and loss of equilibrium might be related, i.e., that magnetic nonequilibrium 
might result from loss of equilibrium, but this must be left as an open 
question. 

The topic of true interest is how an arcade evolves subsequent to 
reaching the critical point. The most appropriate model in this case becomes 
a dynamical one, in which the time-dependent behavior is followed 
explicitly. The dynamical approach resolves many questions that have arisen 
from previous studies of sequences of equilibria. In Section 5.2 we illustrate 
how magnetic nonequilibrium arises during the dynamical evolution of an 
arcade in response to photospheric flow. 

The main emphasis of this paper is to present a detailed study of the 
arcade disruption within the ideal MHD model and to indicate its 
relationship to previous theoretical work. However, we do discuss the 
general effect of finite plasma resistivity on the disruption in order to contrast 
it with ideal MHD behavior, and to emphasize the necessity of finite 
resistivity for impulsive energy release. The resistive behavior will be 
investigated more thoroughly in a future paper. In Sections 2, 3, and 4, we 
describe the plasma model, including the initial and boundary conditions and 
the plasma parameters. In Section 5 we describe our results, and in Section 6 
we discuss the implications for the solar corona. 

2. MODEL DESCRIPTION 

Our motivation is to extend previous studies of the dynamical 
evolution of coronal arcades (Mikic et al. 1988; Biskamp & Welter 1989; Finn 
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et al. 1992). We adopt a spherical coordinate system, finite plasma pressure, 
mass evolution, and gravity, and we use nonuniform grid cells to minimize 
the effect of artificial computational boundaries. We solve the following 
time-dependent MHD equations numerically in spherical coordinates: 

VxB=^J, (1) 


„ „ 1 3B 

V X E = , 

C if 


E + “vxB = r; J , 


§£ 

dt 


+ V-(pv) = 0 , 


dt 




= 7 J x B “ + Pg + V (vpVv) , 


ft+V- 

dt 


(pv) = - (y - l)pV-v + H, 


( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 


where B is the magnetic Held intensity, J is the electric current density, E is the 
electric field, v is the plasma velocity, p and p are the plasma pressure and 
mass density, g = -goR^x/r 1 is the gravitational acceleration, g 0 = 2.74 x 10 4 
cm s -2 is the surface gravity, R 0 = 6.96 x 10 5 km is the solar radius, i] is the 
plasma resistivity, v is the kinematic plasma viscosity, and y is the ratio of 
specific heats. In practice, we use the vector potential A, where B = V x A, to 
implement the algorithm. This model is applicable to both ideal MHD, in 
which we set T) = 0, and to resistive MHD, when i) is finite. The heating term 
in the energy equation (6) accounts for heating due to resistive and viscous 
dissipation, H = (y- l)(rjj 2 + vpVv:Vv). Since numerical resolution 
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constraints (see the Appendix) limit us to using enhanced values for the 
resistivity and viscosity compared to coronal values, the associated heating 
rate may be unphysically large. We therefore frequently neglect this heating 
by setting H = 0; this was done for the simulation described in Section 5.1. We 
have verified that the magnetic field evolution is insensitive to the presence 
of this heating term. The plasma beta, defined by (3 = 8 Jtp/B 2 , measures the 
ratio of plasma to magnetic pressure. These equations are solved in 
axisymmetric spherical coordinates (r,0), in the domain { R 0 < r ^ 

0 < 9 <, 1 80°}, where r is the distance from the center of the Sun, and 6 is the 
solar latitude (with 0 = 0 and 6 = 180° representing the North and South 
poles). In this two-dimensional approximation all quantities are assumed to 
be independent of the solar longitude <p. In the future we plan to extend this 
model to problems in which quantities may vary in all three dimensions 
(r,d,<p). The Alfven speed is given by v A = B/^4xp. Distances are normalized 
by R 0 , and the time scale is normalized by the Alfv6n time, t a = R Q /v A , 
corresponding to a length scale R 0 , evaluated at the equator at the base of the 
corona (r = R 0/ d = 90°). 

The principal results in this paper were calculated using these full 
MHD equations. In order to understand the nature of the disruption, and to 
compare directly with previous theory of force-free fields, we have also used a 
subset of this model: in the zero-beta model we neglect the plasma pressure 
and gravitational forces by setting p = 0 and g = 0, and we assume a fixed 
plasma density profile pit). This model approximates the behavior of the 
low-beta coronal plasma, and contains the same magnetic phenomena that 
are present in the full MHD equations. In particular, equilibria in the zero- 
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beta ideal MHD model are force-free (j x B = 0), so that it is possible to 
compare our results directly to previous analytic and numeric .aiculations 
of force-free fields (see Section 5.2). 

3. INITIAL CONDITIONS AND PLASMA PARAMETERS 


The initial state is a plasma in hydrostatic equilibrium with a dipole 
magnetic field. The pressure and density profiles are chosen to be spherically 
symmetric, p = p(r) and p = p(r), with base values at r = R Q given by p Q and p Q/ 
with a uniform entropy, pp~ y = p 0 Po~ y - The dipole magnetic field, 

„ 2B 0 K 0 3 cos 6 n B o Ko 3 sin0 „ „ 

°r - ^3 / °Q - / of = 0 , 

is a potential field (J = 0) outside the Sun, r £ R 0 . This initial state satisfies the 
static equations (1)— (6) with v = 0 and E = 0. The explicit solution for the 
pressure and density is 


P = Po 



(y - UXqRqPo y, _ «o'|l 1/<r * ” 

rpo l r J. 


/ 



The "photospheric" boundary at r = R 0 is to be regarded as the base of the 
corona (just above the transition region), since we do not attempt to resolve 
the photospheric scale height; this assumption is appropriate because we are 
interested in modeling the large-scale (~ R 0 ) behavior of the corona. The 
following boundary values (at r = R 0 ) were used: magnetic field strength at 
the equator B 0 = 2.2 G; mass density p 0 = 1.67 x 10 -16 g cm -3 (corresponding to 
an ion and electron number density n Q = 10 8 cm” 3 ); temperature 
T 0 = 1.4 x 10 6 K (defined by the ideal gas law p Q = 2nJcT 0 )- The ratio of specific 
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heats is chosen as y = 1.05 to produce a temperature which has little variation 
in the corona (Parker 1963). This hydrostatic atmosphere has plasma 
properties similar to the large scale corona, but does not include the solar 
wind. [The solar wind creates a helmet streamer configuration (Pneuman & 
Kopp 1971); the effects of shear motions on helmet streamer configurations 
C e.g ., Linker, Van Hoven, & Schnack 1990; Linker, Van Hoven, & McComas 
1992) are discussed briefly in Section 6, and will be studied in detail in a future 
paper.] Figure 1 shows the variation of the thermodynamic variables in the 
hydrostatic equilibrium. At the equator, at r = R 0/ these parameters give a 
sound speed c s = 157 km s - \ an Alfven speed v A = 470 km s -1 , and (3 = 0.2. At 
the poles, v A = 940 km s -1 and (3 = 0.05. The Alfven time is t a = 1,449 sec (24.2 
minutes). 

The Lundquist number S is defined as the ratio of the resistive 
diffusion time r R = AkR 0 2 / t\c 2 , to the Alfven time t a . For the Spitzer value of 
the plasma resistivity at a temperature of 1.4 x 10 6 K, we find that S = 5 x 10 14 . 
This large value of S in the corona makes the ideal (i.e., zero resistivity) MHD 
behavior of coronal fields a topic of interest. As will be shown by the results, 
our simulations with r\ = 0 approximate ideal behavior very closely. During 
the resistive computations, such a small value of the plasma resistivity 
cannot be treated accurately in a numerical simulation. For the mesh 
resolutions used in this paper we have been able to accurately simulate the 
resistive plasma evolution at the (much-enhanced) uniform resistivity 
corresponding to S = 10 4 . The Appendix contains a discussion of the 
numerical diffusion in the algorithm and the extent to which ideal MHD 
behavior can be modeled by a numerical solution of the MHD equations. For 
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a uniform plasma viscosity, the viscous diffusion time is defined by 
t v = R 0 2 / v. Typically, we use a viscosity corresponding to t v = 10 2 t a , although 
we have used viscosities as small as t v = 10 3 t a . 

4. GEOMETRY AND BOUNDARY CONDITIONS 

In spherical geometry, two boundaries appear in the simulation: the 
physical inner radial boundary at r = R Q and an artificial outer radial boundary 
at r = Rp which we usually place in the range 50-200R o . In order to maintain 
high spatial resolution in particular regions while minimizing the effect of 
the artificial computational boundary we employ a nonuniform mesh in 
which the computational cells vary in size. The cell size varies very 
gradually, so that neighboring cell dimensions are stretched by only a few 
percent per cell. A typical computational mesh is shown in Figure 2. This 
mesh has 200 x 200 r-6 grid points; the outer radial boundary is at r = 50 R Q ; Ad 
varies from 0.27° at the equator to 2.7° at the poles; A r = 0.01R o at the solar 
surface, increasing gradually to Ar ~ 0.03R o at r ~ 2 R 0 , and to A r - 3R 0 at 
r = 50R o - (This mesh was used for the calculations described in Section 5.1.) 
As will be apparent from the results, a mesh with very small cells in certain 
regions is required to adequately resolve the steep gradients that form. While 
such high-resolution grids are desirable for accuracy, they result in stringent 
time step limitations if a traditional explicit time-integration technique is 
used. To overcome this difficulty we employ a semi-implicit time-integration 
scheme in which the time step is limited by accuracy, rather than by 
numerical stability considerations (see the Appendix). 

We now discuss the boundary conditions at the radial boundaries. 
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Physically, the MHD equations (with finite viscosity) require boundary 
conditions on all components of the flow velocity, the tangential electric field, 
and the normal magnetic field. The build-up of magnetic energy in the 
corona is believed to be driven by flows in the photosphere. Because of the 
high photospheric mass density and conductivity, convective motions of the 
plasma in the photosphere tend to move the footpoints of the coronal 
magnetic field. This situation has been modeled with a "line-tied" boundary 
condition (Raadu 1972; Einaudi & Van Hoven 1981; Priest 1982). This line- 
tied boundary condition is implemented at the lower radial boundary by 
using the "frozen-in" condition given by the ideal Ohm's law. 


„ VpxB 

E +-‘ 7 — -0, 


where v p is the specified photospheric boundary velocity. Only the tangential 
components of E are advanced using this ideal Ohm's law at the boundary, 
since the normal component, E r , is advanced self-consistently by the normal 
component of Ohm's law, equation (3), at the boundary (Mikicef al. 1988). 
Since we do not presently include the solar wind, the normal velocity v r is set 
to zero, we set vq = 0, v<p = v# 0 is a specified shear profile, and p and p evolve 
self-consistently. Equation (7) implies that the tangential electric field at the 
boundary is given by Eg = -v^B r /c and E^ = 0, guaranteeing that the normal 
magnetic field (and hence the photospheric flux) is fixed, dB r /dt = 0. Note that 
it is not necessary (nor is it appropriate) to specify boundary conditions on the 
tangential magnetic fields; they evolve self-consistently according to the 
specified flow. This implementation of boundary conditions is facilitated by 
the use of staggered grids. When the resistivity is zero, this boundary 
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condition causes the magnetic field footpoints to move with the applied flow 
velocity (as verified in Section 5.1); with finite resistivity, there is a small 
relative slip between the applied flow velocity and the motion of the 
magnetic field footpoints that is proportional to rj (Miki6 et al. 1988). 

The upper radial boundary conditions are implemented to minimize 
the effect of the boundary on the solution. When solving the full MHD 
equations, we specify boundary conditions based on the gas characteristics 
(e.g., Courant & Friedrichs 1948) with a non-reflecting condi >n for outgoing 
waves (Heastrom 1979). Neglecting the magnetic field in ; . characteristic 
boundary conditions is permissible since (5 is very high there (for the initial 
dipole field, /J ~ 10 5 at r = 50 Rq). The gas characteristics give boundary values 
for the pressure, density, and the normal (radial) velocity. The tangential 
electric field is determined using equation (7), and the tangential velocity is 
adverted using the normal velocity. This boundary condition has been found 
to work well under both subsonic and supersonic outflow conditions — in 
simulations of helmet streamers an ejected plasmoid propagates out r. f the 
upper radial boundary without reflection. In the case of the zero-beta i/IHD 
equations, we use a simple outward advection scheme, in which all quantities 
are adverted at the local radial velocity when the flow is directed outward. 
This ad-hoc condition is an improvement over a solid perfectly conducting 
wall which allows plasma to leave the domain. Since the upper radial 
boundary is placed far from the Sun (r = 200R o is almost at the position of the 
Earth!), we find that the evolution is not sensitive to the exact boundary 
conditions used (see Section 5.3 for a discussion of the quantitative effect of 
the upper radial boundary). In general, it ought to be possible to implement 
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boundary conditions based on the full MHD characteristics (Hu & Wu 1984). 

The internal boundaries at 0 = 0 and 180° are treated using appropriate 
geometrical conditions required of analytic functions. Even though all the 
results presented in this paper have been found to be symmetric with respect 
to the equator 0 = 90°, all the calculations employed a full domain in the 
region 0<, d<, 180°, so that anti-symmetric modes were not explicitly excluded. 

A shearing motion is introduced by specifying the following 
longitudinal flow profile vf at the photosphere: 

vf = t> 0 (f)©exp[(l - © 4 )/4] / (8) 

where 0= (0 - 9O)/A0 m , and v 0 (t) is a function which is used to specify the 
time profile; typically we turn the flow on and off smoothly using a linear 
ramp over 20-50 t a . The width of the profile is controlled by A 6m) the flow is 
anti-symmetric about the equator 0 = 90°, with a minimum v^° = -v 0 at 
0 = 90 - A0 m , and a maximum v^° = v Q at 0 = 90 + A0 m . We have used 
A0 m = 20°, which concentrates the flow near the neutral line (at the equator), 
and localizes the shear within the region - 50°-130°. This flow profile is 
shown in Figure 3, and is practically identical to that used by Steinolfson 
(1991), except that it has continuous derivatives for all 0. For v 0 > 0, the field 
lines near the equator in the Northern hemisphere are moved in the 
-<t> direction, while those in the Southern hemisphere are moved in the 
+<J> direction. Figure 4 shows a three-dimensional view of the sheared field 
lines. This flow is not intended to represent observed large-scale flows; it is 
merely an idealized flow which produces shear near the neutral line. It is 
convenient to quantify the applied shear using the distance 
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ASmax - J Q V 0 (t ) dt' , (9) 

which is seen to be the maximum distance that a field line footpoint moves 
in the longitudinal direction. For the shear profile given by equation (8), the 
field line with footpoints at 6 = 90 ± A0 m therefore has a total displacement 
between its footpoints equal to 2As max along the Sun's surface in the 
<j> direction. 

5. RESULTS 

In this section we describe the evolution of coronal arcades. In 
Section 5.1 we contrast the ideal and resistive MHD evolution, and we show 
that resistive evolution can lead to the disruption of an arcade. In Section 5.2 
we use the zero-beta MHD equations to identify the cause of the disruption, 
and in Section 5.3 we analyze the results using the virial theorem. 

5.1. Arcade Evolution and Disruption 

We first discuss the evolution of an arcade for a typical calculation. 
The photospheric flow velocity v 0 (t) was increased from zero (linearly in 
time) to 1.9 km s -1 from t = 0 to t = 20t a ; it was kept constant until t = 500t a , 
when it was decreased to zero (linearly in time) at t = 520 t a ; the flow 
remained zero for t £ 520t a . The full MHD equations (l)-(6) were solved, 
starting from the initial conditions described in Section 3. The outer radial 
boundary was placed at r = 50Ro. We describe two calculations which contrast 
the difference between the ideal and resistive evolution. Case 1 is an ideal 
MHD run (with tj = 0), which was run until t = 900t a . Case 2 is identical to 
case 1 until t = 520 t a , but a uniform resistivity corresponding to S = 10 4 is 



introduced at t - 520t a (when the photospheric flow has been turned off). 

Figure 5 shows contours of the flux function y = rsin&A# during the 
ideal evolution (case 1). Contours of y correspond to projections of the 
magnetic field lines in the r-Q plane. The field lines rise in response to the 
applied photospheric shear (Figs. 5a-c) as electric current is induced in the 
corona and a longitudinal field develops. As the shear approaches 
As max = 1.8B 0 (Fig. 5 d), the field lines begin to rise dramatically when 
additional shear is applied. This expansion of the upper field lines is 
accompanied by a squeezing of the lower field lines near the equator in the 
region r = 2-3 R 0 , as seen in Figure 5, and a corresponding concentration of the 
azimuthal current density J^. It will be shown in Section 5.2 that this 
sensitivity is indicative of a fundamental change in the properties of the 
equilibrium magnetic field. The state at t = 520 r A and beyond has 
As m ax = 2.0R o - Once the shear is timed off, the configuration does not change 
significantly in the strong-field regions (Figs. 5e-f). 

When a finite resistivity is introduced at t - 520t a (case 2), the 
subsequent resistive evolution is markedly different from the ideal 
evolution. Figure 6 shows the flux contours during the resistive evolution, 
demonstrating that the state with As max = 2.0K o is susceptible to rapid 
magnetic reconnection. The reconnection causes an X-point to form, creating 
an island (O-point) in the flux. The field lines in this "plasmoid" lose their 
connection to the photosphere, causing the plasmoid to travel upward 
(Figs. 6 e-f). (In our two-dimensional model this plasmoid is an axisymmetric 
torus of helical field lines.) The reconnection is accompanied by a current 
sheet and fast flows, resulting in the dissipation of a substantial fraction of the 
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magnetic energy. Contours of the current density /$ are plotted in Figure 7 at 
t = 549.3 f A (corresponding to the flux plot shown in Fig. 6d), showing the 
equatorial current sheet at the reconnection site. Figure 8 shows the poloidal 
( r-0 ) component of the flow in the vicinity of the reconnection site, 
superimposed on the flux contours, at the same instant. The flow is seen to 
jet out of the reconnection region in the familiar manner, with a maximum 
speed that approaches 600 km s -1 . Figure 9 shows the magnetic energy W for 
cases 1 and 2, normalized by the potential field energy, Wpot = B 0 2 R 0 3 /3. Note 
the stark difference between the ideal and resistive results: whereas the ideal 
MHD evolution shows that the magnetic energy remains approximately 
constant after the shear is turned off, the resistive MHD evolution shows that 
a significant fraction of the free magnetic energy is dissipated. Figure 10 
shows the large increase in kinetic energy K during the reconnection. 

We have computed the energy conservation as a check on the accuracy 
of the numerical method. For the ideal MHD run, the total energy (magnetic, 
thermal, and kinetic) is conserved to within -0.2% (when account is taken of 
energy that leaves or enters the domain, dissipation, and work done against 
gravity). There is a large reservoir of thermal energy that does not change 
substantially during the arcade evolution; when measured in terms of Wpot/ 
the energy error is -2.7%. This demonstrates the good energy conservation 
properties of the code, since the magnetic energy increases by 64% of W po t 
during the ideal evolution. Energy conservation deteriorates during the 
resistive rim, as expected, since the resolution is not as good during the 
impulsive phase of the run when fast flows and a current sheet appear; the 
total energy error is -0.7% ( i.e ., -9% of W po t)* These negative energy errors 
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indicate that there is additional (numerical) diffusion in the algorithm, as 
expected from the upwind treatment of advective terms (see the Appendix). 
This numerical diffusion indicates that the effective Lundquist number 
during the impulsive phase is closer to S = 7 x 10 3 , rather than the assumed 
value S = 10 4 . During the resistive run, the magnetic energy that is lost (34% 
of Wpot) and the thermal energy that is lost or flows out of the domain (7%) is 
distributed as follows: 18% is dissipated resistively, 4% is dissipated viscously, 
0.5% goes into kinetic energy, 9.5% does work against gravity, and 9% is 
presumably dissipated by numerical diffusion. [In this run we neglected the 
heating due to resistive and viscous dissipation by setting H = 0 in 
equation (6). Therefore, the dissipated energy does not appear as heat, but is 
lost from the system, though it is accounted for in the energy conservation 
check.] 

As one test of how closely the simulation approximates ideal behavior 
with 77 = 0, and to check the accuracy with which the line-tied boundary 
conditions have been implemented, we have explicitly traced the field lines at 
the end of the ideal MHD evolution at t = 900 r A . Figure 11 shows a 
comparison between the actual longitudinal field line footpoint positions and 
the positions expected due to the applied flow. The applied shear is 
Asmax s 2.0 R 0 , which implies that the maximum longitudinal displacement 
between field line footpoints ought to equal 4.0R o (for the field line with 
footpoints at $ = 70° and 110°), as shown in Figure 11. Note the excellent 
agreement between the expected and actual footpoint positions (maximum 
error < 0.05K o )/ indicating that the line-tied boundary conditions have been 
implemented accurately. When we perform a simulation in which a 
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nonzero 77 is present during the shearing phase, some slippage of the field 
lines does occur (Miki6 et al. 1988), and the magnetic energy build-up is 
slightly smaller. When 77 is nonzero, there is some slippage which is 
proportional to 77. (The slippage velocity is on the order of lO -4 ^ at S = 10 4 , 
which is small compared to the maximum drive velocity.) However, the 
qualitative evolution of the arcade (and, in particular, the disruption) is 
similar to the case shown. 

Whe he applied shear is smaller than the critical shear (which we 
estimate as approximately Asmax ~ i K 0 )/ the olution of the arcade is quasi- 
static; i.e., if at any point we turn off the shear velocity and allow the system 
to relax to an equilibrium, the configuration changes very little. When Z7 0 is 
well below v A , the magnitude of the shear velocity determines only how 
quickly a particular state is reached, but does not affect the nature of the 
equilibrium for a given applied shear. We have verified that when 
As max £ 1.6R 0 , the configuration settles to an equilibrium when 77 = 0. In the 
presence of finite resistivity, the currents in the configuration diffuse at a rate 
that is directly proportional to 77. There is no impulsive dynamical behavior. 

At larger values of the shear, the solutions undergo a dramatic 
transition. The configuration becomes very sensitive to additional shear, and 
the higher field lines rise significantly. At this point the configuration is 
susceptible to magnetic reconnection when resistivity is present, as described 
above. When the applied shear is close to the critical shear the resistive 
disruption proceeds slowly; however, when the critical shear is exceeded 
significantly, the onset of reconnection is more rapid, and the reconnection 
itself becomes increasingly more vigorous. A detailed study of this 
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reconnection (in particular, its rate as a function of S) will be deferred to a 
future investigation. 

These results indicate that, in the presence of resistivity, an isolated 
arcade disrupts when sufficient shear is applied, with significant liberation of 
magnetic energy, resulting in fast flows and ejection of a plasmoid. In the 
next section we investigate the underlying cause of the disruption within the 
context of the ideal, zero-beta MHD model, and we explore its relation to 
previous work on force-free equilibria. 

5.2. Zero-Beta Model 

Solution of the zero-beta resistive MHD equations shows that the 
evolution of the magnetic field is qualitatively similar to that found using the 
full MHD equations, as described in the previous section, indicating that the 
phenomenon leading to the disruption of the arcade is magnetic in nature. 
The evolution for small shear is quasi-static, and the field disrupts beyond a 
critical shear threshold. 

We have generated a sequence of ideal MHD force-free equilibria by 
applying a velocity v 0 = 0.94 km s -1 for different durations to achieve shears 
of As max = 1-2, 1.4, 1.6, 1.8, 2.0, and 2.2R 0 . After the flow velocity was turned 
off, the simulations were run until the field approached a steady state. The 
outer radial boundary was placed at r = 200R o . Since our primary intention in 
this section is to determine force-free equilibria, we chose a (fixed) density 
profile p = p 0 (R 0 /r) 4 , which produces an Alfven speed that is approximately 
uniform for an open field. This density profile is similar in the lower corona 
to that for hydrostatic equilibrium, but falls off more rapidly at larger r, so that 
the equilibration rate of the weak magnetic fields in the outer corona is 
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enhanced. (By using different density profiles we found that the qualitative 
behavior to be described does not depend on the detailed rorm the dens' ty 
profile, although the time scales may be affected by the Alfv£n speea./ 
Figure 12 shows contours of y/ for these equilibria. A transition in the nature 
of the equilibria is evident: when As max ^ 1.8R 0 , all the field lines are dosed, 
but at higher shears "open" field lines appear. [Although these field lines 
may not be strictly open, i.e., they may close at infinity or at a large distance 
from the Sun, they can be considered to be open for all intents and purposes 
(Sturrock 1991).] The corresponding contours of (Fig. 13) show that this 
transition is accompanied by the formation of a current sheet, i.e., a tangential 
discontinuity in the magnetic field; Figure 14 shows that B r approaches a 
function which has a discontinuity on the equator. Note that the magnitude 
of the maximum current density shown in Fig. 13 increases dramatically as 
the field opens; the current density would be infinite at a true discontinuity in 
B r . It is not usually possible to find ideal MHD numerical solutions with 
discontinuities unless special care is taken. In the absence of resistivity, a 
current sheet usually collapses to the mesh size and generates unphysical 
oscillations, terminating the simulation. Fortunately, we have been able to 
exploit the symmetry in the present configuration (the current sheet appears 
exactly on the equator) to carefully formulate the solution across the equator, 
allowing B r to approach a discontinuous solution (see the Appendix). 

It is important to note that we have not proved conclusively that the 
magnetic field becomes discontinuous. In our numerical scheme we cannot 
rule out the possibility that the magnetic field changes over a narrow layer (of 
finite, rather than zero, width), although it is not readily apparent what 
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physical process would set the width of this layer. The fact that the field 
approaches the open field (which has a true discontinuity in the magnetic 
field) leads us to adopt the hypothesis that the field approaches a true 
discontinuity. 

The "open-field" configuration has been introduced as the asymptotic 
state of highly sheared force-free equilibria (Barnes & Sturrock 1972; Aly 1984, 
1991; Yang et al. 1986; Sturrock 1991). For a given force-free field, the 
corresponding open field is defined as the field with the same photospheric 
flux, but in which the ends of all field lines that intersect the photosphere 
extend to infinity [see Barnes & Sturrock (1972) for a specific example]. The 
magnetic field in the open configuration is potential everywhere except in 
isolated regions in which current sheets (tangential discontinuities in the 
magnetic field) separate oppositely directed field lines. The energy in an 
open-field configuration is an upper limit to the magnetic energy for all force- 
free equilibria (Aly 1984, 1991; Sturrock 1991). For the particular case of a 
dipole flux distribution, we have computed the open field using the 
expedient (involving a monopole field) suggested by Barnes & Sturrock 
(1972). This open field has a discontinuity in B r at the equator, and has a 
magnetic energy W 0 p en = 1.662W po t. The field lines in the open field are 
compared to those in the equilibrium with As ma x = 2.2R 0 in Figure 15. 

Figure 16 shows the height of several field lines in the force-free 
equilibria as a function of applied shear. The field line height increases 
steadily with shear, but becomes very sensitive to shear when a critical shear 
is approached, corresponding to the transition described above. Several 
aspects of this transition should be noted. First, although a fundamental 
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change occurs in the magnetic configurations for the higher shear cases, there 
is no impulsive release of magnetic energy when resistivity is absent. Beyond 
a critical shear, in the ideal MHD model the field lines slowly expand 
outward, opening the field and creating a current sheet (Figs. 12 and 13). 
Second, the reconnection which occurs during the resistive MHD evolution 
is directly related to the transition in the ideal equilibria. There is no 
disruption below the critical shear, whether or not resistivity is present. 
Third, although a current sheet forms in the cases with As max = 2.0 and 2.2 R 0 , 
the magnetic energy remains below the open-field energy throughout the 
evolution. Figure 17 shows the magnetic energy in the force-free equilibria as 
a function of applied shear. Note that magnetic energy approaches the open- 
field energy for large shear, and that the field begins to open (at As max = 2.0 
and 2.2 R 0 ) while W is still below W ope n/ so that these fields represent partially 
open magnetic configurations. The case with As max = 2.2 R 0 has closed field 
lines corresponding to ~ 3% of the flux (these are the low-lying field lines 
near the equator in Fig. 15a). We have made conservative error estimates of 
W for the cases with As max = 2.0 and 2.2R 0 (indicated by the error bars in 
Fig. 17), since the fields have not reached a complete steady state by the end of 
the runs. Our results for the magnetic energies of force-free equilibria are 
therefore in agreement with those of Aly (1984, 1991) and Sturrock (1991). 

It is illuminating to examine the evolution of the toroidal (B^) and 
poloidal (B r , Bq) components of the magnetic field in response to the applied 
shear. Figure 18 shows the evolution of the "partial magnetic energies" W r , 
Wg, and W$, where we have defined W r = Js r 2 dV /Bn, and similarly for Wg 
and W^, for the case As max = 2.0 R 0 . Initially, as the field is sheared, poloidal 
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field is converted into toroidal field (i.e., Wq decreases and W^ increases); 
however, as additional shear is applied, the expansion of the arcade causes 
W# to decrease (Aly 1990), and W r to increase, as the field becomes more 
radial. The partially open field has a small residual toroidal field (in the 
closed-field region). For comparison, note that the fully open field has 
W r = 1.597Wpot, = 0.065 Wpot, and W* = 0. 

The presence of gravity in the full MHD equations would be expected 
to impede the opening of the field. This is evident from Figures 5 and 12, 
which show that whereas a shear of As m ax = 2.0R o is sufficient to open the 
field in the zero-beta model, the corresponding full MHD configuration (with 
gravity) has closed field lines. When the shear is large, the full (ideal) MHD 
model shows that the configuration also approaches an open field. (This is of 
purely theoretical interest, since it is the resistive evolution that is relevant to 
the Sun, and we have already shown in Section 5.1 that the arcade disrupts 
within the full MHD model.) 

5.3. Application of the Virial Theorem 

The transition in behavior past the critical shear can be analyzed using 
the scalar virial theorem (e.g., Aly 1984). A force-free field in equilibrium 
must satisfy the following equality: 

* g2 r / p2 \ 

~ 2 dV = -J I (Bn)(Br) - -y (rn)ldS , (10) 

where V is the volume of the domain, S is the surface bounding the volume, 
and n is the unit outward normal. We refer to the right-hand-side of 
equation (10) as the surface term. In our spherical domain with boundaries at 
r = R 0 and r = Rj we can thus write: 
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f B 2 J 

W ■ — dV = S(R 0 ) - S(R\) , (11) 

Jv 8 tt 

where 

r 3 f/r 

s( r) = j j*(B r 2 - B fl 2 - B* 2 ) sin0d0 . (12) 

How closely W and the surface term match each other is a measure of how 
close a configuration is to equilibrium at a given time. Figure 19 shows the 
evolution of W and the surface term for the case with As max = 0R o . During 
the early part jf th j evolution, W and the surface term match closely, 
verifying that the evolution is quasi-static, but at t ~ 925 r A (when Asmax begins 
to exceed 1.8R 0 ), the two curves start to diverge, indicating that the system is 
no longer in equilibrium. This has been identified as the onset of magnetic 
nonequilibrium (Section 5.2). The two curves approach each other as the 
system equilibrates, although the equilibration time is long. Figure 20, which 
shows the evolution of the height of selected field lines at the equator, 
indicates that the upper field lines rise slowly as the field opens. The long 
equilibration time can be attributed to the slow expansion of the weak upper 
magnetic fields in a region with small Alfv6n speed. 

The expansion which occurs during the opening of the field in the full 
MHD model is also characterized by a long equilibration time, as a result of 
the small Alfv6n speed in the outer corona. It takes thousands of Alfv6n 
times for the field to equilibrate (1000 t a corresponds to 16.8 days), implying 
that the quasi-static approximation is no longer valid once the critical shear is 
reached, even for small photospheric flow speeds. (The field equilibrates on 
the local Alfv6n time scale, which can be much longer in the weak-field 
regions in the outer corona than the "characteristic Alfv4n time scale" t a .) 
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An accurate estimate of the equilibration time depends on the 
thermodynamic model of the outer corona. In particular, a solar wind 
solution has a smaller density in the outer corona than that predicted by the 
hydrostatic model. An important consequence of this slow equilibration is 
that nonequilibrium magnetic fields with energies exceeding W open may 
form in the corona. We have verified that W ope n can be exceeded in the full 
MHD model when the arcade is driven continuously with a photospheric 
shear velocity v 0 = 0.94 km s _1 . 

In order to illustrate the approach to equilibrium when the critical 
shear is not exceeded. Figure 21 shows W and the surface term for the case 
As max = 1.8 R 0 . The evolution is seen to be quasi-static at all times, with the 
field settling to equilibrium relatively quickly when the shear flow is turned 
off at t = 950 t a , as shown by the evolution of the field line height (Fig. 22). 
The residual difference (~5x lO^Wpot) that remains between the W and 
surface term curves in Figures 19 and 21 is a measure of the differencing error 
involved in the computation of the various terms. 

The quantitative influence of the upper radial boundary on the 
solution can be measured by comparing the relative sizes of the three terms 
in equation (11). For the case As max = 2.0R O/ we find that S(R\) remains 
negligibly small during the quasi-static evolution, increasing as the field 
begins to open. Even then, its maximum value (at t ~ 1700t a ) is only 
lO^W^t (i.e., W and S(R 0 ) balance to within 1 part in 10 3 ), indicating that it 
plays a minimal role in the energetics of the field. This demonstrates that the 
upper radial boundary is sufficiently far away so as to have a negligible effect 
on the computed solution. 



28 


6. DISCUSSION 

The results presented in detail in the preceding section can be 
summarized as follows. A sheared arcade responds to slow photospheric 
flows quasi-statically, evolving through sequences of equilibria with 
increasing magnetic energy. With increasing shear, the field lines rise 
dramatically. When a critical shear is exceeded the arcade suffers magnetic 
nonequilibrium: the arcade expands outward from the Sun in a process that 
opens magnetic field lines and forms a current sheet. The current sheet 
results from a tangential discontinuity in the magnetic field which separates 
the two regions of opposite polarity. In the ideal MHD model this opening of 
the field is not very impulsive, and does not produce a large change in 
magnetic energy. It is characterized by a long equilibration time, since it 
involves expansion of weak outer magnetic fields. In marked contrast, when 
finite resistivity is present there is a violent disruption of the configuration 
resulting from magnetic reconnection at the current sheet, leading to 
Alfv6nic flows, considerable dissipation of magnetic energy, and the ejection 
of a plasmoid. This disruption occurs in both the full MHD model and the 
zero-beta model; therefore, it is a property of the magnetic configuration of 
the arcade. 

For the force-free equilibria found in the zero-beta model, the magnetic 
energy remains below the open-field energy (Aly 1984, 1991; Sturrock 1991). 
The critical shear at which magnetic nonequilibrium occurs corresponds to a 
magnetic energy W cr it ~ 1.62W po t, which is slightly (but definitely) below the 
open-field energy, W 0 pen " 1.662Wp 0 t* Thus, following magnetic 
nonequilibrium, the field reaches a partially open configuration in which the 
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majority of the field lines open, with a small fraction of the low-lying field 
lines remaining closed. When the critical shear is exceeded the arcade 
evolution is not quasi-static, indicating that a sequence of equilibria may not 
model coronal evolution accurately during this phase. 

In a study of ideal MHD force-free equilibria of an axisymmetric model 
of a sunspot, Barnes & Sturrock (1972) and Yang et al. (1986) did not identify 
the critical behavior described here. As we have noted, within the ideal MHD 
model the consequences of magnetic nonequilibrium are not very dramatic — 
the field expands slowly toward an open configuration, with the eventual 
formation of a current sheet The equilibration time for this expansion can be 
long compared to the characteristic Alfv6n time at the base of the arcade. [The 
critical behavior may be even less apparent in a numerical calculation of the 
equilibrium field: the current sheet may be artificially broadened by the finite 
resolution of the grid; the expansion of the field may be inhibited by the 
proximity of computational boundaries; and the long equilibration time may 
result in slow convergence toward the open field. These complications were 
recognized by Porter et al. (1992): the large uncertainties in the magnetic 
energy estimates of computed force-free equilibria were attributed to limited 
spatial resolution and the proximity of computational boundaries.] It is 
entirely possible that this sunspot field may also experience magnetic 
nonequilibrium when the magnetic field energy approaches W 0 pe n . It may 
prove worthwhile to revisit the study of the ideal and especially the resistive 
MHD properties of this configuration in light of the present results. 

The disruptive behavior exhibited by a sheared arcade is not unlike 
that seen in previous computations of arcade evolution in Cartesian 
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coordinates. In a very idealized model (which had an artificial array of 
periodic arcades), it was found that arcades disrupted when a critical shear was 
exceeded (Mikic et al. 1988). In that case the mechanism which initiated the 
disruption was ideal MHD instability in a highly symmetric equilibrium, 
although a similar effect was also observed in nonsymmetric arcades 
(Biskamp & Welter 1989). Finn et al. (1992) have discussed the effect of 
symmetry on the relationship between instability and loss of equilibrium in 
coronal arcades. These studies have shown that coronal arcades can disrupt 
through resistive reconnection of magnetic field lines in the presence of 
plasma resistivity. Biskamp & Welter (1989) concluded that "single" arcades 
did not experience disruptive behavior in Cartesian geometry ( i.e ., with 
translational symmetry), in contrast with the disruption observed in groups 
of adjacent arcades, although this conclusion was not based on an exhaustive 
parameter study. In contrast, we have shown that a "single" arcade can 
indeed disrupt. In the (infinitely periodic) configuration studied by Mikic e t 
al. (1988), the corresponding open-field energy is infinite (even when 
measured per arcade, per unit length in the ignorable direction); this is a 
direct consequence of the assumption of translational symmetry. It is thus 
not possible to interpret that disruption as the approach of the magnetic 
energy towards the open-field limit. 

A significant disagreement between our results and those of 
Steinolfson (1991) should be noted. Since our dynamical model is similar to 
that used by Steinolfson, our results ought to give a comparable description of 
arcade evolution. The differences in the results must be attributed to 
differences in the numerical solution of the same physical equations. 
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Steinolfson notes that the evolution of an arcade is very different depending 
on which particular scheme is used to advance the tangential field at the 
photospheric boundary (in one case the field disrupts and in another it does 
not). As noted in Section 4, a boundary condition on B^ is not allowed in a 
proper formulation; when line-tying is specified properly, B^ must be 
evolved self-consistently at the boundary. We have verified that our line- 
tied boundary condition has been implemented accurately (see Figure 11). 
The nature of the eruptive behavior reported by Steinolfson is different from 
the disruption we have described: Steinolfson (1991, Figs. 3, 8, and 9) shows 
that the magnetic energy increases by many orders of magnitude when the 
field disrupts, whereas we have shown that the magnetic field is the source of 
energy for the disruption, so that the magnetic energy decreases (as shown in 
Fig. 9). Steinolfson does not identify the nature of his disruption, nor does he 
explain the source of energy; in contrast, we have explicitly traced the cause of 
the disruption to magnetic nonequilibrium (Section 5.2), and we have 
checked that energy is conserved (Section 5.1). Therefore, the physical 
relevance of Steinolfson's result must be questioned. 

The analytic MHD properties of an idealized model of a coronal arcade 
with an embedded filament have been studied previously (Priest & Forbes 
1990; Forbes & Isenberg 1991). Forbes and Isenberg (1991) found that a 
filament can experience loss of equilibrium in ideal MHD when the flux in 
the filament is varied, although the energy release associated with this event 
is small (less than ~ 1% change in magnetic energy). Forbes & Isenberg (1991) 
and Forbes (1991) have recognized that it is not necessarily the ideal MHD 
energy release associated with loss of equilibrium that drives a coronal mass 
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ejection; it is the subsequent reconnection that releases a substantial amount 
of magnetic energy. As we have demonstrated, once a current sheet forms 
even a small plasma resistivity is effective in dissipating magnetic energy and 
causing plasmoid ejection. Finn & Guzdar (1993) have examined the 
relationship between loss of equilibrium and reconnection in laboratory and 
solar plasmas. 

Although the present model is highly idealized, it is just this type of 
magnet : energy conversion and plasmoid ejection that is required tc '.plain 
coronal mass ejections. In order to make a detailed comparison with 
observations of coronal mass ejections it will be necessary to further refine 
this model. In particular, we will add the important effect of the solar wind 
(Parker 1963). The solar wind opens some of the outer field lines, creating a 
helmet streamer configuration (Pneuman & Kopp 1971; Linker et al. 1990, 
1992). Our preliminary results indicate that a similar disruption also occurs 
in a sheared helmet streamer configuration. In this regard, we conjecture that 
coronal mass ejections are most frequently observed to occur in helmet 
streame: onfigurations because the field is alre iy partially open in such a 

configuration (by the effect of the solar wind), so that a smaller amount of 
shear may be required to disrupt the field. In future work we will study in 
detail the effect of shear on helmet streamers. 

It does not appear to be necessary for the coronal arcade to have 
dimensions comparable to the solar radius for this disruption to occur (as is 
the case in the model problem described here). We therefore speculate that 
active-region-sized arcades ought to exhibit similar behavior. 

In closing, we note that this example offers a good illustration of the 
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mechanism of topological dissipation through magnetic nonequilibrium 
(Parker 1972). It shows that continuous long-wavelength displacements of 
field line footpoints of a magnetic field configuration (which initially has no 
neutral points) eventually lead to the formation of a tangential discontinuity, 
as a consequence of constraints imposed by ideal MHD equilibrium. In the 
presence of finite resistivity the resulting current sheet is seen to be the site of 
rapid reconnection, and leads to dissipation of magnetic energy. 
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APPENDIX 

In this appendix we describe some details of the numerical simulation 
code. The MHD equations (1)— (6) are integrated using finite differences in 
space and time. A leapfrog time integration scheme is used, combined with a 
semi-implicit method. The semi-implicit method allows the time step to 
exceed the Courant limit, which can make the computation more efficient 
than one employing an explicit method. In the semi-implicit method 
(Harned and Kemer 1985; Schnack et al. 1987; Mikic et al. 1988) the time step is 
chosen according to accuracy constraints, rather than by numerical stability 
limitations. Unconditional stability (to wave-like modes) is achieved by 
introducing an implicit linear operator in the momentum equation which 
modifies the dispersion properties of high-Jt modes (where k is the spatial 
wave vector). These are the modes which are not represented accurately by 
the spatial differencing. This scheme is particularly efficient when the slow 
evolution involves long wavelengths (which are treated accurately), but a 
high-resolution mesh is present to capture the strong gradients that form 
during a subsequent impulsive phase. This is precisely the situation that 
occurs in the arcade evolution. The field evolves quasi-statically in response 
to slow, long-wavelength photospheric flows, but eventually a current sheet 
forms, and, in the resistive case, reconnection produces large flows. A 
relatively large time step can thus be used during the quasi-static phase, with 
a corresponding reduction during the impulsive phase. The (explicit) 
advection time step limit provides a natural control of the time step, since it 
automatically reduces the time step when flows become large. The size of the 
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semi-implicit term depends on the size of the time step; it vanishes when the 
time step satisfies the Courant condition. In addition, since the semi-implicit 
term involves a time derivative, it does not affect the steady-state solutions of 
the equations. The matrices arising from the semi-implicit and fully implicit 
viscosity operators are inverted using a preconditioned conjugate gradient 
method (Mikic and Morse 1985). 

As a particular example, for the mesh used in Section 5.1 (shown in 
Fig. 2), the explicit Cour~nt limit for Alfv6n waves is At < 0.004t a . Thus, to 
integrate the quations . 2000t a (as was necessary) would require 500,000 

time steps. Clearly, the evolution does not suggest that this level of accuracy 
is required. With the semi-implicit method we were able to use a time step in 
the range At = 0.025-0.1 r A during the quasi-static evolution, with sufficient 
accuracy and a substantial savings of computer time. We did not find 
significant differences in selected computations when we varied the time 
step. 

Staggered spatial meshes, in which the various fields are defined at 
different (staggered) locations with respect to each other, were employed. 
Such meshes facilitate the implementation of boundary conditions, since it is 
not necessary to specify quantities that do not require a boundary condition, 
and preserve the vector identities V V x A = 0 and V x = 0 in discretized 
form, leading to strict separation of longitudinal and transverse components 
of a vector. In particular, this implies that V B =0 is maintained to within 
roundoff error during the entire computation. 

We now discuss how numerical diffusion is controlled in the code. 
When numerical solutions of the ideal MHD equations are attempted, the 
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following question inevitably arises: how accurately is it possible to maintain 
the constraints imposed by 77 = 0? The answer depends to a large extent on the 
design of the simulation code and the physics of the problem being studied. 
In general, if the evolution generates gradients that are steeper than can be 
resolved on the mesh, the numerical solution will cease to be valid: typically, 
oscillations will develop, causing the simulation to terminate. Sometimes it 
is preferable to avoid this termination by introducing numerical diffusion to 
prevent the gradients from steepening beyond that allowed by the mesh 
resolution, in which case the solution does not represent the ideal evolution, 
but one with small, but finite, resistivity. 

Numerical diffusion in our code can be introduced by using upwind 
(rather than centered) differences of advective terms. Typically we use 
upwind differences when advecting velocity, pressure, and mass, which 
introduces numerical diffusion that has the effect of diffusing steep gradients 
in momentum, mass, and energy. When studying ideal MHD we use 
centered differences to advect the vector potential, to prevent the 
introduction of numerical diffusion in Ohm's law. As a consequence, the 
ideal MHD constraints prohibiting reconnection are preserved. However, if 
the ideal evolution causes a current sheet to form, the simulation will 
terminate unless a specialized treatment is used in the regions where the 
magnetic field becomes discontinuous. When we use the resistive model, we 
find that it is sometimes beneficial to use upwind differencing when 
advecting the vector potential, which introduces a small amount of 
additional diffusion in Ohm's law. This is the source of most of the energy 
error discussed in Section 5.1. 
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In the present arcade evolution, the current sheet that forms with the 
onset of magnetic nonequilibrium occurs exactly on the equator, and is 
characterized by a discontinuity in B r ; however, the symmetry of the situation 
can be exploited. The field approaches an open field, in which this 
discontinuity is in force balance: although B r changes sign discontinuously 
across the equator (i.e., B r jumps from -B e to +B e ), the magnetic pressure B e 2 
remains continuous across the equator. This implies that the part of the 
0-component of the J x B force which equals -(B r 'r)dB r /dd ought to be 
differenced in the form -(l/2r)5B r 2 / 50, with a choi of staggering in which 
the B r mesh-points straddle (but do not fall on) the equator. When the 
equations are differenced in this form, the approach to a discontinuous 
solution can be accommodated in the ideal MHD model without introducing 
numerical diffusion in Ohm's law. 

That this technique can successfully handle a current sheet on the 
equator is illustrated by the stark difference between the ideal and resistive 
evolution (see Fig. 9), which indicates that the reconnection is controlled by 
the physical resistivity. We have found that in the ideal model, numerically 
induced reconnection is absent when this technique (and sufficient 
resolution) is used. 

The code has been designed to run on supercomputers (we used Cray- 
YMP and Cray-YMP/C90 machines), and takes anywhere from a few minutes 
to many hours of CPU time, depending on the parameters. 
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FIGURE CAPTIONS 

Fig. 1. — Radially symmetric profiles of pressure p, mass density p, and 
temperature T, in the initial hydrostatic equilibrium. The base values at 
r = R 0 are given in the text. 

Fig. 2. — The mesh cells for a typical calculation. The upper radial 
boundary L- at r = 50R o - The i.. it shows the detail of the mesh - the solar 
surface. 

Fig. 3. — The normalized longitudinal shear velocity profile v<p° that is 
applied at the photosphere. 

Fig. 4. — A three-dimensional view of the sheared magnetic field lines for 
Asmax = 1-6R 0 - At f = 0, in the dipole field, the selected field lines lie in the 
(r-6) plane. 

Fig. 5. — Evolution of the projected field lines (contours of constant flux 
Y) in the ideal MHD case. The applied shear at .he times indicated in the 
sequence (a)-(/) is As max = 0, 1.0, 1.4, 1.8, 2.0, and 2.0R o . (The shear flow is 
turned off at t = 520t a .) The field lines in the initial dipole field are shown 
in (a). Note the sensitivity of the field line height to shear as As ma x 
approaches 1.8R 0 - The contour values in this and all subsequent y^-plots 
correspond to 19 uniformly spaced values between 0.05 Yo and 0.95 Yo* 
inclusive, where Yo = B 0 R 0 2 is the maximum flux. 

Fig. 6. — Evolution of the projected field lines during the resistive MHD 
case. The reconnection of field lines creates an O-point in the flux and ejects a 
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plasmoid. 

Fig. 7. — Contours of the longitudinal component of the electric current 
density at t = 549.3 r A , corresponding to the flux contours shown in Fig. 6 d 
for the resistive case. A current sheet appears at the equator, with a 
maximum current density equal to 31.2(cB 0 /4^R 0 ). 

Fig. 8. — Detail of the reconnection region at t = 549.3t a , showing the 
poloidal projection of the flow vectors, with the field lines superimposed. 
The maximum flow speed approaches 600 km s _1 . 

Fig. 9. — Evolution of the magnetic energy W for the ideal and resistive 
MHD cases. The shear flow is turned off at t = 520t a . In the resistive case, 
reconnection dissipates a large fraction of the free magnetic energy. In 
contrast, in the ideal case the magnetic energy changes little after the shear is 
turned off. The solid symbols and open circles indicate the values of time at 
which the flux is shown in Figs. 5 and 6, respectively. 

Fig. 10. — Evolution of the kinetic energy K for the ideal and resistive 
MHD cases. The shear flow is turned off at t = 520 t a . The kinetic energy 
increases considerably during the reconnection. The solid symbols and open 
circles indicate the values of time at which the flux is shown in Figs. 5 and 6, 
respectively. 

Fig. 11. — Comparison of the computed and exact longitudinal field line 
footpoint displacements at the end of the ideal MHD rim at t = 900 r A . The 
solid line shows the exact footpoint displacements expected from the applied 
flow profile; the symbols denote the displacements computed by tracing 
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individual field lines. The agreement shows that the line-tied boundary 
conditions have been implemented accurately. 

Fig. 1Z — Contours of the flux y for ideal force-free MHD equilibria as a 
function of applied shear. The outer field lines open when As max ~ 2.0Ro. 

Fig. 13. — Contours of the longitudinal component of the electric current 
density corresponding to the flux contours shown in Fig. 12 for ideal force- 
free MHD equilibria. The maximum is expressed in units of the normalizing 
current density, cB 0 / 4 jcR 0 . A current sheet appears when Asmax >Z0Ro- 

Fig. 14. — Radial magnetic field B r in the force-free field with 
As max = 2.2R 0 (corresponding to Figs. 12 f and 13 ft at r = 3R 0 / showing that B r 
approaches a function which has a discontinuity at the equator 6 = 90°. This 
is the location of the current sheet in the open field. The symbols show the 
locations of the Q mesh points. 

Fig. 15. — Comparison of the projected field lines (contours of the flux y) 
for (a) the force-free field with shear As ma x = 2.2 R 0 , and (b) the open field 
configuration. The force-free field corresponds to a partially open field, with a 
small fraction of closed low-lying field lines near the equator. The contour 
values are the same in both plots, and correspond to the usual 19 uniformly 
spaced values of y; a contour at y= 0.98y o has been added in both plots to 
emphasize the closed field lines in (a). 

Fig. 16. — Height of selected field lines in the ideal MHD equilibria as a 
function of the applied shear. The field lines are labeled by the corresponding 
flux value y, in terms of the maximum flux y Q . The field line height rises 
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sharply as the critical shear (As max - 1.8R 0 ) is approached. For Asmax > I.ORq, 
as the field opens, the higher field lines (with y = 0.95, 0.9, and 0.8 1^ 0 ) are 
characterized by h -» <». 

Fig. 17. — Magnetic energy W of force-free equilibria as a function of 
applied shear Asmax- Note that W approaches the open-field energy for large 
shear, but that the field opens while W is still below W 0 pen (at Asmax = 2.0 and 
2.2 R 0 ). The error bars indicate uncertainties in W since the indicated states 
are not in a complete steady state. 

Fig. 18. — Evolution of the "partial magnetic energies" W r/ Wg, and W^, 
which are defined as the energies in the r, 0, and <p components of the 
magnetic field, for the case As ma x = 2.0R o . The shear flow is turned off at 
t = 1050t a . 

Fig. 19. — Evolution of the magnetic energy W and the surface term in 
the virial theorem comparison for the case As max = 2.0R o . The two curves 
match closely during the quasi-static evolution. The magnified view in the 
inset shows that the curves depart when magnetic nonequilibrium occurs. 
The subsequent equilibration occurs over a long time scale. The "bump" in 
the surface term at t ~ 1600t a occurs when the expanding upper field lines 
reach the upper radial boundary. 

Fig. 20. — Evolution of the height of selected field lines at the equator 
(0 = 90°) for the case As max = 2.0 R 0 . The field lines are labeled by the 
corresponding flux value y, in terms of the maximum flux y Q . The shear 
flow is ramped down to zero starting at t = 1000t a and ending at t = 1050t a , 
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with zero flow thereafter. When the critical shear is exceeded the upper field 
lines rise to large heights as the field begins to open. The slow asymptotic rise 
of the upper field lines indicates a long equilibration time. The evolution of 
the low field lines indicates that there is a closed-field region in the resulting 
partially open field. 

Fig. 21. — Evolution of the magnetic energy W and the surface term in 
the virial theorem comparison for the case As max = 1.8R 0 . The magnified 
view in the inset shows that the curves remain close at all times, indicating 
quasi-static evolution, in contrast with the case when the critical shear is 
exceeded (Fig. 19). 

Fig. 22. — Evolution of field line height for the case As^x = 1.8R 0 . When 
the shear flow is ramped down to zero, starting at t * 900t a and ending at 
t = 950 t a , the field lines equilibrate quickly, in contrast to the behavior 
observed when the critical shear is exceeded (Fig. 20). 
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